# An automated method for recovering piecewise smooth functions on spheres free from Gibbs oscillations.

Abstract

Spectral methods using spherical harmonic coefficients are widely used in applications for geophysics and atmospheric sciences. The major drawback in using spherical harmonic spectral methods occurs when the underlying function is piecewise smooth. In this case, the well-known Gibbs phenomenon reduces the order of accuracy to first order and produces spurious oscillations, particularly ill regions near the discontinuities. The frequently studied Gegenbauer reconstruction method has been shown to alleviate the effects of the Gibbs phenomenon while restoring the exponential accuracy of the spectral approximation. Since each reconstruction must be implemented only within smooth regions, the jump discontinuities of the piecewise smooth function must first be located by an edge detection method. This study combines the recent developments in both edge detection and Gegenbauer reconstruction methods to construct an automated procedure for recovering horizontal or two-dimensional geophysical global fields free from Gibbs oscillations and without compromising the high order convergence properties of spectral methods. Numerical efficiency and robustness are discussed.

Key words and phrases : Gibbs phenomenon, Gegenbauer polynomials, spherical harmonics, edge detection.

2000 AMS Mathematics Subject Classification--42A10, 42C05, 42C20, 65M70.76E20, and 861A10.

1 Introduction

Spectral methods using the two-dimensional spherical harmonic basis are frequently used for solving problems that occur on spheres [6, 22, 28]. For example, they are commonly employed in geophysics to describe lateral changes on the Earth since the data sets are sparse and also since polynomials make a convenient choice in describing the Earth's functions. They are also used to reconstruct seismic images where the data can be extremely sparse. Another application is in weather forecasting models that solve the atmospheric equation of motion where, in addition to their high convergence properties, spherical harmonic spectral methods have proven to be quite capable of simulating the long range behavior of the climate system [26]. We note that finite difference methods are also popular choices for such applications, and have the advantage of being simpler to formulate computationally. In this paper we discuss how to improve some of the difficulties that arise from using spherical harmonic spectral methods.

Spherical harmonic spectral methods have the same limitations as every other polynomial expansion method. Specifically, if the underlying function has jump discontinuities, then the approximation suffers from artificial oscillations appearing near the jump discontinuities, and the overall convergence rate is reduced to first order. This behavior is known as the Gibbs phenomenon [7, 18]. Unfortunately in the case of geophysical and atmospheric models, the Gibbs phenomenon is quite prevalent because the data are sparsely accumulated and what should appear as steep gradients are actually seen as jump discontinuities in this under-resolved environment. The "Gibbs ringing artifact" is widespread in these locations (see, e.g., Figure 5.2). Such ringing can have disastrous consequences. For instance, in the case of an atmospheric model coupled to an ocean model, the wind's spurious Gibbs oscillations in the mountain field (caused by the lack of resolution of the steep gradients) has a nonlinear interaction with the ocean model. The effects of the oscillations in the model are intensified and over time these oscillations will lead to complete instability of the model. As another example, in the case of reconstructing seismic images, the oscillations can be misidentified as genuine artifacts. Important characteristic features of the true image will be misunderstood causing all results based on the reconstruction to be flawed.

In effort to combat the Gibbs phenomenon for both the modeling and reconstruction of geophysical data, some sort of smoothing filter is often employed [7, 27]. Filtering is a powerful tool in modeling because it restores high accuracy away from the discontinuities while maintaining stability. However, since filtering causes the finer small scale features of a horizontal or two-dimensional geophysical global field to be "smoothed" over, long term simulations can be severely impacted as their solutions can be heavily dissipated in time. Similarly, filtering in the geophysical global field can cause a blurring effect resulting in the loss of important information in the data.

Other techniques, such as the spectral viscosity method [25, 32], have been designed to simulate long-term nonlinear models with as little damping as possible to maintain the integrity of the solution coefficients. Specifically, the method allows some oscillations to occur as long as the overall stability is not jeopardized. Once the simulation is complete the underlying piecewise smooth solution must still be recovered. It was demonstrated in [30] that the spectral viscosity method yields solution coefficients contain enough information for high resolution recovery of the underlying solution when applied to one-dimensional conservation laws. The method was adapted to solve the shallow water equations on a sphere with spherical harmonic basis functions [13], and it is anticipated that analogous results can be shown for the spectral viscosity spherical harmonic solution coefficients.

The purpose of this study is to provide an automated high-order resolution reconstruction method that will recover a horizontal or two-dimensional geophysical global field without Gibbs oscillations or undesirable blurring. The method should be suitable both to post-process spectral viscosity solutions as well as to reconstruct fields from spherical harmonic coefficients, e.g., the mountain fields in weather forecasting problems. We wish to avoid filtering that causes "blurring" near the discontinuities and destroys the finer features of the solution. In this preliminary investigation we consider the reconstruction of piecewise smooth functions from their spherical harmonic coefficients. The immediate benefits can be seen in the reconstruction of the mountain fields for climatology problems. In addition we anticipate that the method will work well to post-process the numerical solutions of partial differential equations on a sphere, provided that the solution retains enough high mode information, i.e., the solution should not be overly damped. Such a solution allows the Gibbs ringing to remain within the simulation while maintaining stability, and the reconstruction is applied only after the final time step. We do not, however, tackle problems that have irregularities near a pole. This will be saved for future investigations.

One popular method for recovering piecewise smooth functions from their spectral coefficients is the Gegenbauer reconstruction method, which was introduced in [21] and adapted for spherical harmonics in [12]. The advantage of the Gegenbauer reconstruction method over other types of reconstruction algorithms is that exponential accuracy is recovered up to the jump discontinuities. This is particularly important in applications in geophysics and climatology, since much of the activity occurs near the jump discontinuities. Other methods have been effectively used to combat the Gibbs phenomenon when the given information is Fourier data [8, 9, 29]. So far, however, none of these methods have been explored for spherical harmonics.

Recall that all high-resolution reconstruction methods, including the Gegenbauer reconstruction method, must avoid crossing jump discontinuities and therefore can only be used in smooth regions. Hence, it is necessary to first determine where the jump discontinuities of the piecewise smooth functions are, as their placement defines the boundary points of each smooth region. An edge detection technique that computes jump discontinuities from spectral coefficients was developed in [16]. For the sake of simplicity and automation, in this paper we adopt the method designed in [3] which determines jump discontinuities from local physical grid point data. Although the method is designed to determine jump discontinuities for multivariate functions, for our purposes here it is only necessary to apply the method dimension by dimension. A side benefit is that no complicated two-dimensional algorithm needs be employed. The Gegenbauer reconstruction method can then be used to restore a horizontal or two-dimensional geophysical global field in each smooth sub-region, and as a final step these sub-region reconstructions are "glued together" to restore--Gibbs-free--either the original data field or the underlying spectral viscosity solution from a long term simulaton.

There have been many studies conducted on the Gegenbauer reconstruction and edge detection methods. The purpose of this study is to unify the approaches and to construct an automated procedure that is robust, easily implemented, and cost-efficient for horizontal or two-dimensional geophysical global fields. This paper is organized as follows: The problem of reconstructing a piecewise smooth function f([theta], [PHI]) on a sphere is formulated in [section] 2. The local edge detection algorithm, which determines the regions of smoothness for the function, is described in [section] 3. In [section] 4 the Gegenbauer reconstruction method is first reviewed for one-dimensional problems using both Fourier and orthogonal polynomial expansions. It is then suitably modified for the spherical ease. Issues of numerical implementation are also discussed. Examples in [section] 5 are presented to confirm that the spherical harmonic coefficients contain enough information to obtain a highly accurate reconstruction to a global data field in a robust and flexible way. Finally, in [section] 6 we discuss our results and future applications.

2 Preliminaries

To mathematically formulate the problem, we begin by defining the spherical harmonic spectral expansion of a function f([theta], [PHI]) with colatitude coordinate [theta] [member of] [0, [pi]] and longitude coordinate [PHI] [member of] [0, 2[pi]] as follows:

Definition 2.1.

f ([theta], [phi]) = [[infinity].summation over (q=0) [summation over ([absolute value of v] [less than or equal to] q) [a.sup.v.sub.q][Y.sup.v.sub.q]([theta], [phi]) (21)

where the spherical harmonics [Y.sup.v.sub.q] ([theta], [phi]) are given by

[Y.sup.v.sub.q]([theta], [phi]) = [square root of (2q + 1)(q-v)!/4[pi](q + v)!][P.sup.v.sub.q](cos[theta])[e.sup.iv[phi]].

Recall that [P.sup.v.sub.q] (cos[theta]) are the associated Legendre polynomials orthogonal on [-1, 1]. The coefficients [a.sup.v.sub.q] are given by

[a.sup.v.sub.q] = [[integral].sup.2[pi].sub.0] [[integral].sup.[pi].sub.0]] [integral] ([theta],[phi]) [[Y.sup.v.sub.q]([theta], [phi])]* sin [theta]d[theta]d[phi], (2.2)

where [[Y.sup.v.sub.q]([theta], [phi])]* are the complex conjugates of [Y.sup.v.sub.q]([theta], [phi]).

The standard reconstruction of a function is performed by taking the truncated spectral representation of f([theta], [phi]),

fN([theta], [phi]) = [N.summation over(q-0} [summation over ([absolute value v][less than or equal to] q) [a.sup.v.sub.q][Y.sup.v.sub.q]([theta], [[phi]). (2.3)

As stated in the introduction, if f([theta], [phi]) is smooth over the entire sphere, then f N([theta], [phi]) converges exponentially to f([theta], [phi]). However, if the underlying function has jump discontinuities, then the Gibbs phenomenon will occur. Moreover, this behavior will also occur if the steep gradients of f([theta], [phi]) are not properly resolved, as is often the ease in geophysical and atmospheric models. The following example, displayed in on a rectangular domain [0, [pi] X [0, 2[pi]] in Figure 2.1, shows the undesirable consequences of applying (2.3) to reconstruct a piecewise smooth function.

Example 2.1.

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2.4)

Filtering is often used to alleviate the effects of the Gibbs ringing artifact. One major drawback in using filters, as seen ill Figure 2.1(c), is that the finer features of a function are often lost along with the Gibbs oscillations. It is also possible that some ringing artifacts will remain. Hence, we seek to construct a high order reconstruction algorithm that rids the function of Gibbs oscillations while maintaining the high order convergence properties of the spherical harmonic partial stun expansion for smooth functions.

[FIGURE 2.1 OMITTED]

3 The Local Edge Detection Method

Edge detection is a critical first step in all high-resolution image reconstruction methods. As noted previously, the Gibbs phenomenon results from crossing over a jump discontinuity in the domain. In order to reduce the Gibbs effects without the undesirable blurring caused by filtering, it is first necessary to determine each region of smoothness, which are separated by the "edges" of the function.

Although many edge detection algorithms can be utilized for this purpose, we employ the multivariate method introduced in [3] for its simplicity and robustness. As will be explained in [section] 4, the Gegenbauer reconstruction method is performed dimension by dimension in the latitudinal and longitudinal directions. Hence, the local edge detection method used here needs only be applied as a one-dimensional algorithm. In this case, the method is reduced to computing various orders of Newton divided differences and is therefore robust, efficient, and easy to implement.

Assume that we are given the spherical harmonic coefficients (2.2) so that we can compute the spherical harmonic partial sum [f.sub.N]([theta], [phi]) from (2.3) for any ([[theta].sub.i], [[phi].sub.j]) in the domain [0, [pi]] x [0, 2[pi]]. For simplicity, we assume that these points are chosen apriori as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

We wish to determine whether or not there are "edges" at each ([[theta].sub.i], [[phi].sub.j]) on the given domain. Note that the grid points do not have to be equally spaced. To construct the one-dimensional edge detection algorithm, let us first assume that we want to determine the jump discontinuities along a fixed longitudinal coordinate [[phi].sub.f] = [[phi].sub.j]. Then (2.3) is the one-dimensional function approximation [f.sub.N]([theta], [[phi]sub.f]) given on the data set So. To ease notation, let us say that we are seeking to find the jump discontinuities for a given function f([theta]). Clearly, [f.sub.N]([theta]) is not a good approximation of f([theta]), but we use it since it exhibits jump discontinuities at the same grid point values as f([theta]). It is possible that the Gibbs ringing in [f.sub.N]([theta]) might falsely portray additional jump discontinuities. However, the numerical experiments conducted here indicate that these spurious responses do not affect the robustness of the edge detection method.

We begin describing the edge detection method by defining a jump function

[f]([theta]) := f([[theta].sup.+]) - f([[theta].sup.-]), (3.1)

where f([[theta].sup.+]) and f([[theta].sup.-]) are the right and left side limits of the function f at [theta]. If f is continuous at [theta], then the jump function [f]([theta]) = 0. Also, if

[[THETA].sub.*] = {[[theta].sup.*] : 0 [less than or equal to] [[theta].sup.*] < [pi]}

denotes the set of jump discontinuities of f in [0, [pi]] for fixed [[phi].sub.f], then for any [[theta].sup.*] [member of] [[THETA].sup.*] we have [f]([[theta].sup.*])=f([[theta].sup.*+]) - f([[theta].sup.*-]) [not equal to] 0.

For the mth order local edge detection method, we use a local set of m + 1 points around the point [theta],

[S.sub.m] = {[[theta].sub.1], [[theta].sub.2], ..., [[theta].sub.m+1]}. (3.2)

Note that these points do not correspond to the same point values in [S.sub.[[theta], but rather they are the points that surround the particular value [theta] for which we want to determine whether or not a jump discontinuity exists. The local edge detection method is defined as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (3.3)

which is constructed to approximate the jump function [f]([theta]). The coefficients [c.sub.i]([theta]) in (3.3) are designed to annihilate polynomials of m - 1 degree and are found by solving

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3.4)

where [delta] is the Dirac delta function

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

The normalization factor [q.sub.m]([theta]) is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3.5)

and by design [q.sub.m]([theta]) [not equal to] 0 [3]. The construction of (3.3) leads to

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3.6)

where [I.sub.[theta]] is the smallest closed interval such that [S.sub.m] [subset] [I.sub.[theta]], with [S.sub.m] defined in (3.2). The proof of (3.6) is given in [3]. Loosely speaking, (3.6) states that if there are no jump discontinuities occurring in the set [S.sub.m], then in regions of smoothness (away from the jumps) we have [L.sub.m]f([theta]) [right arrow] 0 with mth order accuracy. At the jump discontinuities (3.3) converges to the value of the jump.

As described in [3], the one-dimensional mth order edge detection method can be simplified to building the mth degree Newton's divided difference. Recall that for the local set defined in (3.2), the mth degree Newton divided difference for [S.sub.m] when m [greater than or equal to] 1 is defined by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3.7)

where f[[theta]] [equivalent to] f([theta]). We can compute (3.7) as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3.8)

for

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3.9)

As shown in [3], the edge detection method (3.3) can be expressed using Newton divided differences as

[L.sub.m]f([theta]) = [m!/[q.sub.m]([theta])] f[[S.sub.m]]. (3.10)

Furthermore, an explicit formula for the coefficients [c.sub.i]([theta]) can be derived as

[c.sub.i]([theta]) = m!/[[omega].sub.i]([S.sub.m]), i = 1, ... , m + 1. (3.11)

yielding the values [q.sub.m]([theta]) in (3.5). Details are found in [3].

Consider the following example displayed in Figure 3.1 sampled on 64 Gaussian points. Since the data are sparsely spaced, it is difficult to determine the steep gradients from the true jump discontinuities.

[FIGURE 3.1 OMITTED]

Example 3.1.

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3.12)

The results of the edge detection method are exhibited in Figure 3.2, where it is clear that [L.sub.m]f([theta]) [right arrow] 0 away from the discontinuities, and as predicted, the convergence increases with the order m.

However the application of (3.3) also encounters some problems in the approximation of jump functions. Specifically ms m increases, oscillations that occur in the neighborhood of a jump discontinuity can be misidentified as true edges, as is evident in Figure 3.2(d). On the other hand, for smaller m there is a risk of identifying a steep gradient as an edge, as seen in Figure 3.2(a). These problems can be exacerbated when the true function f([theta]) is not known and the reconstruction [f.sub.N]([theta]) must be used instead.

To avoid the possibility of misidentification of an edge, we adopt the minmod algorithm originally described in [17] and adapted in [3] for the local edge detection method. Jump discontinuities are distinguished from neighborhood oscillations by effectively incorporating information provided by the results from (3.3) using the various orders of m. This technique is referred to in [3] as the minmod local edge detection method.

Definition 3.1. For a given finite set M [subset] of N of positive integers, consider the set [L.sub.M] f = {[L.sub.m] f : R [right arrow] R | m [member of] M}. The mimnod function is defined by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (3.13)

Loosely speaking, the minmod function distinguishes artificial oscillations from existing jump discontinuities by recognizing that a true jump will not change sign for arbitrary m, while artificial oscillations will vary in sign depending on the order m chosen in (3.3). The application of the minmod function

Figure 3.3: The minmod algorithm (3.13) applied to Example 3.1 using M = 5. The original function is sampled on 64 Gaussian spaced grid points.

[FIGURE 3.2 OMITTED]

[FIGURE 3.3 OMITTED]

(3.13) is exhibited in Figure 3.3.

Even after applying (3.13), it is possible that some small oscillations may still persist in the jump function approximation. This is easily rectified by applying a critical threshold parameter that is scaled according to the particular image being analyzed as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3.14)

Here [L.sub.M]f defines the final jump approximation of f([theta]).

We now review the total procedure of the minmod local edge detection method. First, (3.3) is applied for various m. Then (3.13) and (3.14) are used to ensure that the jumps are correctly identified. The edges are then defined as each value [[theta].sub.i] such that [L.sub.M]f([[theta].sub.i]) [not equal to] 0 for i = 1, ... , [N.sub.[theta]]. The threshold value, [T.sub.crit], is problem dependent. However, [T.sub.crit] [approximately equal to] [DELTA][theta] is a reasonable choice assuming that the image is sufficiently resolved for the purpose of determining jump discontinuities. Specifically, we assume that only one jump discontinuity can exist between two neighboring points. If an error is made such that more edges are found than actually exist, then a high order reconstruction method will still work as the function is still smooth in each region of approximation. On the other hand, if an edge goes undetected, then the function will be approximated over a discontinuity resulting in the Gibbs phenomenon. In this case aposteriori error estimates should be used to determine if the edges must be recomputed.

To determine the edges of a piecewise smooth function from spherical harmonic data, we assume that we are given the grid point values f([[theta].sub.i], [[phi].sub.j]) for i = 1, ... , [N.sub.[theta]] and j = 1, ... , [N.sub.[phi]] over the domain [0, [pi]] x [0, 2[pi]]. One dimension is held fixed, e.g., [phi] = [[phi].sub.f], and then the edges are determined for each f([theta], [[phi].sub.f]) as the set

[[THETA].sup.*] = {[[theta].sup.*] : 0 [less than or equal to] [[theta].sup.*] [less than or equal to] [pi]}.

Each edge [[theta].sup.*] then constitutes an end point for a smooth interval of f([theta], [[phi].sub.f]). Specifically if

[[THETA].sup.*] = {[[theta].sup.+], [[theta].sup.*.sub.2], ..., [[theta].sup.*.sub.j]},

then

(0, [[theta].sup.*.sub.1]), ([[theta].sup.*.sub.1], [[theta].sup.*.sub.2]), ([[theta].sup.*.sub.2], [[theta].sup.*.sub.3]), ..., ([[theta].sup.*.sub.j], [pi])

form the smooth subintervals in which f([theta],[[phi].sub.f]) is to be reconstructed. The method is performed in exactly the same way in the longitudinal direction by holding [theta] = [[theta].sub.f] fixed. In this case, the edge detection algorithm is further simplified since the longitudinal data points [[phi].sub.j], j = 0, ..., [N.sub.[phi]] are uniformly distributed [3]. Hence, for any point ([theta], [phi]) in the domain [0, [pi]] x [0, 2[pi]], we can determine the smooth region which encloses it [[[theta].sub.a], [[theta].sub.b]] x [[[phi].sub.a], [[phi].sub.b]]. These will be the intervals used for the Gegenbauer reconstruction described in [section]4.

4 The Gegenbauer Reconstruction Method

In this section we review the basic concepts of the Gegenbauer reconstruction method [19, 20, 21]. Although we will ultimately be considering the reconstruction of a piecewise smooth function from spherical harmonics, it is useful to review the one-dimensional Gegenbauer reconstruction method for both trigonometric and general orthogonal polynomials expansions. Since the spherical harmonic partial stun expansion contains both Fourier and associated Legendre polynomials, the Gegenbauer reconstruction method must be amenable to both types of expansions.

As discussed in [section]1, the Gegenbauer reconstruction method can only be performed in regions of smoothness. Crossing over a discontinuity in a piecewise smooth function will lead to the Gibbs phenomenon. The endpoints of each smooth sub-region for both the latitudinal and longitudinal directions are first determined by the local edge detection method explained in [section]3.

4.1 Review of the Gegenbauer Polynomials

Before we describe the Gegenbauer reconstruction method, it is useful to briefly review the Gegenbauer polynomials and some related properties. More details can be found in [1, 5, 7, 18].

The Gegenbauer polynomials [C.sup.[lambda].sub.l]([xi]) on the interval [-1, 1] are orthogonal under weight function [(1 - [[xi].sup.2]).sup.[[lambda] - 1/2] for [lambda] [greater than or equal to] 0 such that

[[integral].sup.1.sub.1][(1 - [[xi].sup.2]).sup.[lambda]- 1/2][C.sup.[lambda].sub.k]([xi])[C.sup.[lambda].sub.l]([xi])d[xi] = [[delta].sub.l,k][h.sup.[lambda].sub.l],

where

[h.sup.[lambda].sub.l] = [[pi].sup.1/2][C.sup.[lambda].sub.l](1) [GAMMA]([lambda] + 1/2)/[GAMMA]([lambda])(l + [lambda]) and [C.sup.[lambda].sub.l](1) = [GAMMA](l + 2[lambda])/l![GAMMA](2[lambda]).

A function f([xi]) on [-1, 1] can be represented as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

As with other orthogonal polynomial expansions, the truncated approximation

[f.sub.m]([xi]) = [m.summation over (l=0)] [[??].sup.[lambda].sub.l][C.sup.[lambda].sub.l]([xi])

converges exponentially to f([xi]) as long as f([xi]) is smooth [7, 18]. Note that [lambda] - 1/2 and [lambda] = 0 yield the standard Legendre and Chebyshev polynomial expansions, respectively.

To reconstruct f(x) on a sub-interval x [member of] [a, b] we apply the simple linear transformation

x = [epsilon][xi] + [delta] (4.1)

where [epsilon] = b-a/2 and [delta] = b+a/2. In this case the Gegenbauer coefficients,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.2)

yield the approximation

[f.sub.m](x([xi])) = [m.summation over (l=0)][[??].sup.[lambda].sub.l,[epsilon]][C.sup.[lambda].sub.l]([xi]), (4.3)

which exponentially converges to f(x) in the interval [a, b] as m [right arrow] [infinity] when f(x) is smooth. In what follows, we review how the Gegenbauer polynomials can be applied to reconstruct piecewise analytic functions when the given data. prohibits the direct use of (4.2) and (4.3).

4.2 Development of the Gegenbauer Reconstruction Method

Assume that a function f([xi]) is piecewise smooth on an interval [-1,1] and analytic in the region a [less than or equal to] x [less than or equal to] b, where -1 [less than or equal to] a < b [less than or equal to] 1. We wish to reconstruct f(x) in the sub-interval [a, b] given either the Fourier coefficients,

[[??].sub.k] = 1/2 [[integral].sup.1.sub.-1] f([xi])[e.sup.-ik[pi][xi]] d[xi], k [member of] [- N/2, N/2] (4.4)

or orthogonal polynomial coefficients (e.g., Chebyshev or Legendre),

[a.sub.l] = 1/[gamma]l] [[integral].sup.1.sub.-1] f([xi])[P.sub.l]([xi])w([xi])d[xi], l [member of] [0, N], (4.5)

where [P.sub.l]([xi]) are orthogonal polynomials on [-1, 1] with weight function w([xi]) and normalization factor [[gamma].sub.l].

In the Fourier case, the approximation to f(x([xi])) using the transformation (4.1) is given by the partial sum

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.6)

Analogously,

[f.sub.N] (x([xi])) = [N.summation over (l=0)] [a.sub.l][P.sub.l](x([xi])) (4.7)

represents the standard orthogonal polynomial approximation. It is well known that if f([xi]) is smooth (and periodic in the Fourier expansion case) on [-1, 1], then the approximations (4.6) and (4.7) converge spectrally to f(x([xi])) for any interval [a, b]. However, if f is only piecewise smooth, the above approximations yield very poor results with Gibbs oscillations visible near the endpoints a and b as well as an overall reduced convergence throughout the interval [7, 18]. This is due to the fact that, even though the function may be smooth in the interval of reconstruction [a, b], the coefficients (4.4) and (4.5) are obtained by integrating a piecewise smooth function across the entire region [-1, 1], including the discontinuities. Furthermore, the Fourier approximation poses an additional problem since there is no reason to assume that f is periodic in [a, b].

The Gegenbauer reconstruction method, originally developed in [21], addresses this issue by observing that it is possible to reconstruct a function f(x) by reprojecting the approximations (4.6) and (4.7) onto the orthogonal Gegenbauer polynomials to construct (4.3) in a smooth sub-interval [a, b]. Clearly, the coefficients (4.2) are not known for this reconstruction; rather, (4.6) and (4.7) must be used to construct the approximate coefficients

[[??].sup.[lambda].sub.l,[epsilon]] = 1/[h.sup.[lambda].sub.l] [[integral].sup.1.sub.-1] [f.sub.N](x([xi]))[C.sup.[lambda].sub.l]([xi])[(1 - [[xi].sup.2]).sup.[lambda] - 1/2] d[xi] [approximately equal to] [[??].sup.[lambda].sub.l,[epsilon]], l [member of] [0, m]. (4.8)

Note that if (4.8) is derived from the Fourier sum (4.6), then there is an explicit formula [5] given by:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.9)

where [J.sub.l+[lambda]]([epsilon]k[pi]) is the Bessel function of the first kind [1]. If instead (4.8) is constructed from (4.7), we use the Chebyshev Gauss Labotto quadrature rule,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.10)

where the quadrature points are defined as ([[xi].sub.j] = -cos [[pi]j/N. For the coefficients (4.8) we have

q(x([[xi].sub.j])) = [f.sub.N]([epsilon][[xi].sub.j] + [delta])[(1 - [[xi].sup.2.sub.j]).sup.[lambda][C.sup.[lambda].sub.l]([[xi].sub.j])

and [??] [greater than or equal to] 2[lambda]+m+N. Another useful identity is the closed form for the Christoffel Darboux sum, which will be practical in numerically implementing (4.12) if the coefficients (4.8) are obtained from the orthogonal polynomial coefficients (4.5). That is:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.11)

where

[k.sub.m] = [2.sup.m][GAMMA]([lambda] + m)/m![GAMMA]([lambda]).

Once (4.8) is computed, the Gegenbauer reconstruction method

[g.sub.m](x([xi])) = [m.summation over (l=0)] [[??].sup.[lambda].sub.l,[epsilon]][C.sup.[lambda].sub.l]([xi]) (4.12)

approximates (4.3). As shown in [19, 20], (4.12) is an exponentially convergent approximation to f(x) in [a, b] provided A and m grow with N. (Note that this is why neither the Cheybshev or Legendre expansions in (4.7) can be used for the reprojection (4.12).) In this way, the Gibbs phenomenon can be completely removed up to the endpoints of each smooth sub-interval. Hence, the Gegenbauer method is much more powerful than standard filtering methods, since no blurring will occur at the "edges" of each region. In the investigations of [12, 15] the method has been shown to be effective for two-dimensional and spherical problems, as well as for various applications.

Remarks.

1. The method also applies when the given information are the pseudospectral coefficients rather than (4.4) or (4.5). This is particularly useful when the underlying image is given on grid-point data, as is the case for atmospheric spectral transform models.

2. Although theoretically m and [lambda] grow with N, in practice it has been shown that much smaller values will yield high resolution reconstructions. Smaller values of m and [lambda] are also necessary to avoid round off error and to reduce computational cost, [10, 11].

3. The Gegenbauer reconstruction method works well when the information is given as the partial sum expansion of general orthogonal polynomials (4.7), including associated Legendre polynomials, making it amenable to spherical harmonics, [12, 15].

We now wish to adapt the Gegenbauer reconstruction method to reconstruct piecewise smooth functions on spheres. This will be accomplished by reconstructing the grid-point values f([[theta].sub.i], [[phi].sub.j]), i = 1, ..., [N.sub.[theta]] and j = 1, ... , [N.sub.[phi]], over the whole domain [0, [pi]] x [0, 2[pi]]. Typically, the number of resolved grid-points, [N.sub.[theta]] and [N.sub.[phi]], are related to the number of given spherical harmonic coefficients (2.2). To demonstrate the Gegenbauer reconstruction method, we look a single point ([theta], [phi]) for which we want to reconstruct the function f. Let us assume that f([theta], [phi]) is analytic in the region [[theta].sub.a] [less than or equal to] [theta] [less than or equal to] [[theta].sub.b] and [[[phi].sub.a] [less than or equal to] [phi] [less than or equal to] [[phi].sub.b] which has been previously determined in [section]3.

To apply the Gegenbauer reconstruction method from spherical harmonics, we first make the linear transformations

[theta]([[xi].sub.[theta]]) = [[epsilon].sub.[theta]][[xi].sub.[theta]] + [[delta].sub.[theta]] [phi]([[xi].sub.[phi]]) = [[epsilon].sub.[phi]][[xi].sub.[phi]] + [[delta].sub.[phi]] (4.13)

where [[epsilon].sub.[theta]] = [[theta].sub.b]-[[theta].sub.a]/2, [[epsilon].sub.[phi]] = [[theta].sub.] - [[theta].sub.a]/2, [[delta].sub.[theta]] = [[theta].sub.b] + [[theta].sub.a]/2, and [[delta].sub.[phi]] = [[phi].su.b] + [[phi].sub.a]/2. The spherical harmonic partial sum in the region of reconstruction is

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.14)

As discussed in the introduction and visible in Figure 2.1, the spherical harmonic sum approximation (4.14) yields poor results, with Gibbs ringing particularly prevalent near the edges of the regions. However, noting that the spherical harmonic partial sum is a combination of Fourier and orthogonal polynomial expansions, we can apply the Gegenbauer reconstruction method in the same way as discussed above. Specifically, we use (4.14) to construct the coefficients

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.15)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

This leads to the spectrally convergent approximation [15]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.16)

4.3 Implementation of the Gegenbauer Reconstruction Method

In order for the Gegenbauer reconstruction method to be as effective as a postprocessor for geophysical and atmospherics simulations, it is important that the method be numerically efficient, robust, and easy to implement. These issues were first addressed in [11] and the resulting algorithm is briefly reviewed here.

To compute (4.15) efficiently, we first rewrite

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Applying (4.9) yields

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.17)

The Gauss Labotto quadrature formula is used to approximate

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.18)

where [c.sub.j] is defined in (4.10) for [xi][theta],j = -cos[pi]j/N, and [bar.N] [greater than or equal to] 2[[lambda].sub.1] + [m.sub.1] + N. Employing (4.17) and (4.18) enables us to rewrite (4.16) as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.19)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.20)

The direct application of the Christoffel Darboux formula (4.11) yields an explicit formula for the sum

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.21)

Hence, by defining

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.22)

we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.23)

There are several numerical advantages to computing (4.16) as (4.23). First, the Fourier expansion component of spherical harmonics (2.1) means that the explicit formula (4.17) can be applied instead of costly integration, which allows the use of FFT algorithms for the reconstruction in the longitudinal direction. Second, although we are forced to use a quadrature rule (4.10) to compute the integral in the latitudinal direction, the Christoffel Darboux formula (4.11) allows us to eliminate one of the sums in (4.20) which also reduces the cost of computation. Finally, by "splitting" the reconstruction (4.23) into latitudinal and longitudinal approximations, it is clear that the reconstruction can be performed by either fixing the latitudinal coordinate and computing

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.24)

or alternatively by fixing the longitudinal coordinate and computing

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.25)

Hence, the Gegenbauer reconstruction method can be tailored to a particular application. For example, if the global data field has more discontinuities in the longitudinal direction, it is reasonable to hold [theta] fixed and apply (4.24), and similarly for the latitudinal direction using (4.25). Clearly either of these approximations are more cost efficient than (4.23). We also note that it is possible to apply one of the reconstruction algorithms developed in [8, 9, 29] in the longitudinal direction. This study does not apply such a "splitting" technique to reconstruct a function, but this option will be explored in future investigations.

5 Numerical Examples

To demonstrate the effectiveness of the Gegenbauer reconstruction method for spherical harmonics, we recall Example 2.1. After determining the regions of smoothness by applying the local edge detection method in both the longitudinal and latitudinal directions, the Gegenbauer reconstruction method (4.16) is used to obtain a high-order reconstruction of the image without compromising the finer features of the function. Figure 5.1 shows the results. Note that the parameter choices [[lambda].sub.1] = [[lambda].sub.2] = 4, and [m.sub.1] = [m.sub.2] = 8 are not optimized, but rather chosen for ease of implementation while still ensuring convergence. Parameter optimization and the reduction of possible round off error for the Gegenbauer reconstruction method are discussed in [10, 14, 23]. As an example closer to real applications, we look at a picture of the surface mountains of the Earth. In this ease, although there are no actual discontinuities in the image, the resolution used is insufficient to properly resolve the steep gradients. Hence, as displayed in Figure 5.2(b), the Gibbs oscillations are visible in the model mountains field as a series of altitude oscillations over the oceans. They are especially prominent in the vicinity of sharp mountain ridges like the Andes. As discussed in [section]1, the Gibbs oscillations ill the mountain field is one of the causes of the instabilities arising in the atmospheric model coupled to the ocean model. Although filtering alleviates some of the ringing effects, it is clear from Figure 5.3(a) that the regions near the sharp mountain ridges are blurry, and some oscillations remain. The Gegenbauer reconstruction method (4.16), based on the spherical harmonic coefficients (2.2) for its reconstruction, is exhibited in Figure (5.3), where it is clear that the Gibbs oscillations are removed and accuracy is restored.

[FIGURE 5.1 OMITTED]

[FIGURE 5.2 OMITTED]

Remarks. There are several factors to be considered here.

[FIGURE 5.3 OMITTED]

1. The Gegenbauer reconstruction method should not be used in very small regions that contain too few grid points, as resolution requirements prohibit spectral reconstruction in this case. Local interpolation methods or splines could be used (with diminished accuracy in those regions) and will be discussed in future work. In these examples we used the filtered reconstruction values in the small regions.

2. As discussed in [section]4, it is also possible that one-dimensional reconstruction in either the latitudinal or longitudinal direction would yield sufficient accuracy. This would make the Gegenbauer reconstruction method less costly, particularly if the reconstruction was performed for longitudinal values since the FFT can be applied. Aposteriori error estimates can determine whether or not the full Gegenbauer reconstruction method should be applied, and will be discussed in future work.

3. No attempts at optimization were made, nor was a check of round off error enforced. This has been done in other applications [10, 14, 23], and could be useful here as well.

4. We have not attempted our algorithms for problems that have irregularities near the poles. This can cause added complications, since the reprojection method requires large enough smooth regions (with at least [pi] points per interval) in order to do the reconstruction. A similar issue was addressed in [2]. There the Gegenbauer reconstruction method was used to recover MRI brain images from Fourier data. In small smooth regions, where there were too few points to construct the Gegenbauer reconstruction, a

simple line segment was drawn across the interval. Perhaps something more sophisticated can be used on the sphere when the singularities occur near the pole.

6 Concluding Remarks

This study combines the edge detection algorithm [3], the minmod procedure [17], and the Gegenbauer reconstruction method, e.g., [19, 21, 12], to produce an automated high order reconstruction method for piecewise smooth functions using the spherical harmonics partial sum approximation (2.3). The numerical examples provided in [section] 3 and [section] 5 demonstrate the effectiveness of the method. By incorporating numerical techniques such as Newton divided differencing, FFT, Gaussian quadrature, and the Christoffel Darboux sum, tim method can be made numerically robust and cost efficient. Future work will include the application of the reconstruction method to real geophysical and climatology data, where it will be necessary to develop techniques that validate the reconstructed solutions. In addition, the method will be incorporated within the simulation of models.

References

[1] M. Abramovitz and I. Stegun, Handbook of Mathematical Functions, Dover Publications, New York, 1965.

[2] R. Archibald and A. Gelb, A Method to Reduce the Gibbs Ringing in MRI Scans While Keeping Tissue Boundary Integrity, IEEE Transactions on Medical Imaging, 21:4, 305-319, 2002.

[3] R. Archibald, A. Gelb, and J. Yoon, Polynomial Fitting for Edge Detection in irregularly Sampled Signals and Inmges, SIAM Journal on Numerical Analysis, 43, 259-279, 2005.

[4] R. Archibald, A. Gelb, and J. Yoon, Determining the Locations of the Discontinuities in the Derivatives of Functions, in press, 2006.

[5] H. Bateman, Higher Transcendental Functions, Vol. 2, McGraw-Hill, New York 1953.

[6] J.P. Boyd, Chebyshev and Fourier Spectral Methods, Springer-Verlag, 1989.

[7] C. Canuto, M.Y. Hussaini, A. Quarteroni, and T.A. Zang, Spectral Methods in Fluid Dynamics, Springer 1987.

[8] T.A. Driscoll and B. Fornberg, A Pade-based Algorithm for Overcoming the Gibbs Phenomenon, emphNumerical Algortihms, 26, 77-92, 2001.

[9] K. S. Eckhoff, On a High Order Numerical Method for Functions with Singularities, Mathematics of Computation, 67, 1063-1087, 1995.

[10] A. Gelb, Parameter Optimization and Reduction of Round Off Error for the Gegenbauer Reconstruction Method, Journal of Scientific Computing, 20:3, 433-459, 2004.

[11] A. Gelb, A Hybrid Approach to Spectral Reconstruction of Piecewise Smooth Functions, Journal of Scientific Computing, 15, 293-322, 2001.

[12] A. Gelb, The Resolution of Gibbs Phenomenon for Spherical Harmonics, Mathematics of Computation 66:218, 699-717, 1997.

[13] A. Gelb and J. Gleeson, Spectral Viscosity for Shallow Water Equations in Spherical Geometry, Monthly Weather Review, 129:9, 2346-2360, 2001.

[14] A. Gelb and Z. Jackiewicz, Determining Analyticity for Parameter Optimization of the Gegenbauer Reconstruction Method, SIAM J. Sci. Comput 27:3, 1014-1032, 2006.

[15] A. Gelb and A. Navarra, Recovering Grid-Point Values without Gibbs Oscillations in Two Dimensional Domains on the Spheres, Center for Research on Parallel Computation, California Institute of Technology, 1997.

[16] A. Gelb, and E. Tadmor, Detection of Edges in Spectral Data II. Nonlinear Enhancement, SIAM J. Namer. Anal. 38, 1389-1408, 2000.

[17] A. Gelb, and E. Tadmor, Adaptive Edge Detectors for Piecewise Smooth Data Based on the MinMod Limiter, Journal of Scientific Computing, 28-23, 279-306, 2006.

[18] D. Gottlieb and S. Orszag, Numerical Analysis of Spectral Methods: Theory and Applications. SIAM, Philadelphia, 1977.

[19] D. Gottlieb and C.-W. Shu, A General Theory for the Resolution of the Gibbs Phenomenon, w, Tricomi's Ideas and Contemporary Applied Mathematics Convention, Roma Accademia Nazionale Dei Lincei, Italy, 1998.

[20] D. Gottlieb and C.-W. Shu, On the Gibbs Phenomenon and its Resolution, SIAM Review, 1997.

[21] D. Gottlieb, C.-W. Shu, A. Solomonoff, and H. Vandeven, On the Gibbs Phenomenon I: Recovering Exponential Accuracy from the Fourier Partial Sum of a Non-Periodic Analytic Function, J. Comput. Appl. Math., 43, 81-92, 1992.

[22] G.J. Haltiner and R.T. Williams, Numerical Predictions and Dynamics Meteorology, Wiley, 1980.

[23] Z. Jackiewicz, Determination of optimal parameters for the Chebyshev-Gegenbauer reconstruction method, SIAM J. Sci. Comput., 25:4, 11871198, 2004.

[24] C. Lindberg and A. Broccoli, Representation of topography in spectral climate models and its effect on simulated precipitation, J. Climate, 9, 26412659, 1996.

[25] Y. Maday, S.M. Ould Kaber, and E. Tadmor, Legendre Pseudospectral Viscosity Method for Nonlinear Conservation Laws, SIAM J. Num. Anal., 30, 321-342, 1993.

[26] F.B. Mitchell, The Greenhouse Effect and Climate Change, Rev. of Geophys., 27, 115-139, 1989.

[27] A. Navarra, W.F. Stern, and K. Miyakoda, Reduction of the Gibbs Oscillations in Spectral Model Simulations, J. Climate, 7, 1169-1183, 1994.

[28] S.A. Orszag, Fourier Series on Spheres, Mon. Weather Rev. 102:56, 1974.

[29] B.D. Shizgal and J.H. Jung, Towards the Resolution of the Gibbs Phenomena, Journal of Computational and Applied Mathematics, 161, 41-65, 2003.

[30] C.W. Shu and P. Wong, A note on the accuracy of spectral method applied to nonlinear conservation laws, Journal of Scientific Computing, 10, 357369, (1995).

[31] P.N. Swarztrauber, On the Spectral Approximation of Discrete Scalar and Vector Functions on a Sphere, Siam J. Numer. Anal., 16, 934-949, 1979.

[32] E. Tadmor, Super Viscosity and Spectral Approximations of Nonlinear Conservation Laws, Numerical Methods for Fluid Dynamics IV, M. J. Baines and K. W. Morton eds., Clarendon Press, Oxford, 1993.

Chris Blakely

The University of Maryland, CSCAMM

4146 CSIC Building # 406, Paint Branch Drive, College Park, MD 20742-3289

cblakely@cscamm.umd.edu

Anne Gelb

Department of Mathematics and Statistics, Arizona State University

Tempe, Arizona 85287

ag@math.la.asu.edu

Antonio Navarra

National Institute of Geophysics and Volcanology (INGV)

Via Donato Creti 12. [-40128 Bologna, Italy

navarra@bo.ingv.it

Spectral methods using spherical harmonic coefficients are widely used in applications for geophysics and atmospheric sciences. The major drawback in using spherical harmonic spectral methods occurs when the underlying function is piecewise smooth. In this case, the well-known Gibbs phenomenon reduces the order of accuracy to first order and produces spurious oscillations, particularly ill regions near the discontinuities. The frequently studied Gegenbauer reconstruction method has been shown to alleviate the effects of the Gibbs phenomenon while restoring the exponential accuracy of the spectral approximation. Since each reconstruction must be implemented only within smooth regions, the jump discontinuities of the piecewise smooth function must first be located by an edge detection method. This study combines the recent developments in both edge detection and Gegenbauer reconstruction methods to construct an automated procedure for recovering horizontal or two-dimensional geophysical global fields free from Gibbs oscillations and without compromising the high order convergence properties of spectral methods. Numerical efficiency and robustness are discussed.

Key words and phrases : Gibbs phenomenon, Gegenbauer polynomials, spherical harmonics, edge detection.

2000 AMS Mathematics Subject Classification--42A10, 42C05, 42C20, 65M70.76E20, and 861A10.

1 Introduction

Spectral methods using the two-dimensional spherical harmonic basis are frequently used for solving problems that occur on spheres [6, 22, 28]. For example, they are commonly employed in geophysics to describe lateral changes on the Earth since the data sets are sparse and also since polynomials make a convenient choice in describing the Earth's functions. They are also used to reconstruct seismic images where the data can be extremely sparse. Another application is in weather forecasting models that solve the atmospheric equation of motion where, in addition to their high convergence properties, spherical harmonic spectral methods have proven to be quite capable of simulating the long range behavior of the climate system [26]. We note that finite difference methods are also popular choices for such applications, and have the advantage of being simpler to formulate computationally. In this paper we discuss how to improve some of the difficulties that arise from using spherical harmonic spectral methods.

Spherical harmonic spectral methods have the same limitations as every other polynomial expansion method. Specifically, if the underlying function has jump discontinuities, then the approximation suffers from artificial oscillations appearing near the jump discontinuities, and the overall convergence rate is reduced to first order. This behavior is known as the Gibbs phenomenon [7, 18]. Unfortunately in the case of geophysical and atmospheric models, the Gibbs phenomenon is quite prevalent because the data are sparsely accumulated and what should appear as steep gradients are actually seen as jump discontinuities in this under-resolved environment. The "Gibbs ringing artifact" is widespread in these locations (see, e.g., Figure 5.2). Such ringing can have disastrous consequences. For instance, in the case of an atmospheric model coupled to an ocean model, the wind's spurious Gibbs oscillations in the mountain field (caused by the lack of resolution of the steep gradients) has a nonlinear interaction with the ocean model. The effects of the oscillations in the model are intensified and over time these oscillations will lead to complete instability of the model. As another example, in the case of reconstructing seismic images, the oscillations can be misidentified as genuine artifacts. Important characteristic features of the true image will be misunderstood causing all results based on the reconstruction to be flawed.

In effort to combat the Gibbs phenomenon for both the modeling and reconstruction of geophysical data, some sort of smoothing filter is often employed [7, 27]. Filtering is a powerful tool in modeling because it restores high accuracy away from the discontinuities while maintaining stability. However, since filtering causes the finer small scale features of a horizontal or two-dimensional geophysical global field to be "smoothed" over, long term simulations can be severely impacted as their solutions can be heavily dissipated in time. Similarly, filtering in the geophysical global field can cause a blurring effect resulting in the loss of important information in the data.

Other techniques, such as the spectral viscosity method [25, 32], have been designed to simulate long-term nonlinear models with as little damping as possible to maintain the integrity of the solution coefficients. Specifically, the method allows some oscillations to occur as long as the overall stability is not jeopardized. Once the simulation is complete the underlying piecewise smooth solution must still be recovered. It was demonstrated in [30] that the spectral viscosity method yields solution coefficients contain enough information for high resolution recovery of the underlying solution when applied to one-dimensional conservation laws. The method was adapted to solve the shallow water equations on a sphere with spherical harmonic basis functions [13], and it is anticipated that analogous results can be shown for the spectral viscosity spherical harmonic solution coefficients.

The purpose of this study is to provide an automated high-order resolution reconstruction method that will recover a horizontal or two-dimensional geophysical global field without Gibbs oscillations or undesirable blurring. The method should be suitable both to post-process spectral viscosity solutions as well as to reconstruct fields from spherical harmonic coefficients, e.g., the mountain fields in weather forecasting problems. We wish to avoid filtering that causes "blurring" near the discontinuities and destroys the finer features of the solution. In this preliminary investigation we consider the reconstruction of piecewise smooth functions from their spherical harmonic coefficients. The immediate benefits can be seen in the reconstruction of the mountain fields for climatology problems. In addition we anticipate that the method will work well to post-process the numerical solutions of partial differential equations on a sphere, provided that the solution retains enough high mode information, i.e., the solution should not be overly damped. Such a solution allows the Gibbs ringing to remain within the simulation while maintaining stability, and the reconstruction is applied only after the final time step. We do not, however, tackle problems that have irregularities near a pole. This will be saved for future investigations.

One popular method for recovering piecewise smooth functions from their spectral coefficients is the Gegenbauer reconstruction method, which was introduced in [21] and adapted for spherical harmonics in [12]. The advantage of the Gegenbauer reconstruction method over other types of reconstruction algorithms is that exponential accuracy is recovered up to the jump discontinuities. This is particularly important in applications in geophysics and climatology, since much of the activity occurs near the jump discontinuities. Other methods have been effectively used to combat the Gibbs phenomenon when the given information is Fourier data [8, 9, 29]. So far, however, none of these methods have been explored for spherical harmonics.

Recall that all high-resolution reconstruction methods, including the Gegenbauer reconstruction method, must avoid crossing jump discontinuities and therefore can only be used in smooth regions. Hence, it is necessary to first determine where the jump discontinuities of the piecewise smooth functions are, as their placement defines the boundary points of each smooth region. An edge detection technique that computes jump discontinuities from spectral coefficients was developed in [16]. For the sake of simplicity and automation, in this paper we adopt the method designed in [3] which determines jump discontinuities from local physical grid point data. Although the method is designed to determine jump discontinuities for multivariate functions, for our purposes here it is only necessary to apply the method dimension by dimension. A side benefit is that no complicated two-dimensional algorithm needs be employed. The Gegenbauer reconstruction method can then be used to restore a horizontal or two-dimensional geophysical global field in each smooth sub-region, and as a final step these sub-region reconstructions are "glued together" to restore--Gibbs-free--either the original data field or the underlying spectral viscosity solution from a long term simulaton.

There have been many studies conducted on the Gegenbauer reconstruction and edge detection methods. The purpose of this study is to unify the approaches and to construct an automated procedure that is robust, easily implemented, and cost-efficient for horizontal or two-dimensional geophysical global fields. This paper is organized as follows: The problem of reconstructing a piecewise smooth function f([theta], [PHI]) on a sphere is formulated in [section] 2. The local edge detection algorithm, which determines the regions of smoothness for the function, is described in [section] 3. In [section] 4 the Gegenbauer reconstruction method is first reviewed for one-dimensional problems using both Fourier and orthogonal polynomial expansions. It is then suitably modified for the spherical ease. Issues of numerical implementation are also discussed. Examples in [section] 5 are presented to confirm that the spherical harmonic coefficients contain enough information to obtain a highly accurate reconstruction to a global data field in a robust and flexible way. Finally, in [section] 6 we discuss our results and future applications.

2 Preliminaries

To mathematically formulate the problem, we begin by defining the spherical harmonic spectral expansion of a function f([theta], [PHI]) with colatitude coordinate [theta] [member of] [0, [pi]] and longitude coordinate [PHI] [member of] [0, 2[pi]] as follows:

Definition 2.1.

f ([theta], [phi]) = [[infinity].summation over (q=0) [summation over ([absolute value of v] [less than or equal to] q) [a.sup.v.sub.q][Y.sup.v.sub.q]([theta], [phi]) (21)

where the spherical harmonics [Y.sup.v.sub.q] ([theta], [phi]) are given by

[Y.sup.v.sub.q]([theta], [phi]) = [square root of (2q + 1)(q-v)!/4[pi](q + v)!][P.sup.v.sub.q](cos[theta])[e.sup.iv[phi]].

Recall that [P.sup.v.sub.q] (cos[theta]) are the associated Legendre polynomials orthogonal on [-1, 1]. The coefficients [a.sup.v.sub.q] are given by

[a.sup.v.sub.q] = [[integral].sup.2[pi].sub.0] [[integral].sup.[pi].sub.0]] [integral] ([theta],[phi]) [[Y.sup.v.sub.q]([theta], [phi])]* sin [theta]d[theta]d[phi], (2.2)

where [[Y.sup.v.sub.q]([theta], [phi])]* are the complex conjugates of [Y.sup.v.sub.q]([theta], [phi]).

The standard reconstruction of a function is performed by taking the truncated spectral representation of f([theta], [phi]),

fN([theta], [phi]) = [N.summation over(q-0} [summation over ([absolute value v][less than or equal to] q) [a.sup.v.sub.q][Y.sup.v.sub.q]([theta], [[phi]). (2.3)

As stated in the introduction, if f([theta], [phi]) is smooth over the entire sphere, then f N([theta], [phi]) converges exponentially to f([theta], [phi]). However, if the underlying function has jump discontinuities, then the Gibbs phenomenon will occur. Moreover, this behavior will also occur if the steep gradients of f([theta], [phi]) are not properly resolved, as is often the ease in geophysical and atmospheric models. The following example, displayed in on a rectangular domain [0, [pi] X [0, 2[pi]] in Figure 2.1, shows the undesirable consequences of applying (2.3) to reconstruct a piecewise smooth function.

Example 2.1.

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2.4)

Filtering is often used to alleviate the effects of the Gibbs ringing artifact. One major drawback in using filters, as seen ill Figure 2.1(c), is that the finer features of a function are often lost along with the Gibbs oscillations. It is also possible that some ringing artifacts will remain. Hence, we seek to construct a high order reconstruction algorithm that rids the function of Gibbs oscillations while maintaining the high order convergence properties of the spherical harmonic partial stun expansion for smooth functions.

[FIGURE 2.1 OMITTED]

3 The Local Edge Detection Method

Edge detection is a critical first step in all high-resolution image reconstruction methods. As noted previously, the Gibbs phenomenon results from crossing over a jump discontinuity in the domain. In order to reduce the Gibbs effects without the undesirable blurring caused by filtering, it is first necessary to determine each region of smoothness, which are separated by the "edges" of the function.

Although many edge detection algorithms can be utilized for this purpose, we employ the multivariate method introduced in [3] for its simplicity and robustness. As will be explained in [section] 4, the Gegenbauer reconstruction method is performed dimension by dimension in the latitudinal and longitudinal directions. Hence, the local edge detection method used here needs only be applied as a one-dimensional algorithm. In this case, the method is reduced to computing various orders of Newton divided differences and is therefore robust, efficient, and easy to implement.

Assume that we are given the spherical harmonic coefficients (2.2) so that we can compute the spherical harmonic partial sum [f.sub.N]([theta], [phi]) from (2.3) for any ([[theta].sub.i], [[phi].sub.j]) in the domain [0, [pi]] x [0, 2[pi]]. For simplicity, we assume that these points are chosen apriori as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

We wish to determine whether or not there are "edges" at each ([[theta].sub.i], [[phi].sub.j]) on the given domain. Note that the grid points do not have to be equally spaced. To construct the one-dimensional edge detection algorithm, let us first assume that we want to determine the jump discontinuities along a fixed longitudinal coordinate [[phi].sub.f] = [[phi].sub.j]. Then (2.3) is the one-dimensional function approximation [f.sub.N]([theta], [[phi]sub.f]) given on the data set So. To ease notation, let us say that we are seeking to find the jump discontinuities for a given function f([theta]). Clearly, [f.sub.N]([theta]) is not a good approximation of f([theta]), but we use it since it exhibits jump discontinuities at the same grid point values as f([theta]). It is possible that the Gibbs ringing in [f.sub.N]([theta]) might falsely portray additional jump discontinuities. However, the numerical experiments conducted here indicate that these spurious responses do not affect the robustness of the edge detection method.

We begin describing the edge detection method by defining a jump function

[f]([theta]) := f([[theta].sup.+]) - f([[theta].sup.-]), (3.1)

where f([[theta].sup.+]) and f([[theta].sup.-]) are the right and left side limits of the function f at [theta]. If f is continuous at [theta], then the jump function [f]([theta]) = 0. Also, if

[[THETA].sub.*] = {[[theta].sup.*] : 0 [less than or equal to] [[theta].sup.*] < [pi]}

denotes the set of jump discontinuities of f in [0, [pi]] for fixed [[phi].sub.f], then for any [[theta].sup.*] [member of] [[THETA].sup.*] we have [f]([[theta].sup.*])=f([[theta].sup.*+]) - f([[theta].sup.*-]) [not equal to] 0.

For the mth order local edge detection method, we use a local set of m + 1 points around the point [theta],

[S.sub.m] = {[[theta].sub.1], [[theta].sub.2], ..., [[theta].sub.m+1]}. (3.2)

Note that these points do not correspond to the same point values in [S.sub.[[theta], but rather they are the points that surround the particular value [theta] for which we want to determine whether or not a jump discontinuity exists. The local edge detection method is defined as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (3.3)

which is constructed to approximate the jump function [f]([theta]). The coefficients [c.sub.i]([theta]) in (3.3) are designed to annihilate polynomials of m - 1 degree and are found by solving

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3.4)

where [delta] is the Dirac delta function

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

The normalization factor [q.sub.m]([theta]) is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3.5)

and by design [q.sub.m]([theta]) [not equal to] 0 [3]. The construction of (3.3) leads to

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3.6)

where [I.sub.[theta]] is the smallest closed interval such that [S.sub.m] [subset] [I.sub.[theta]], with [S.sub.m] defined in (3.2). The proof of (3.6) is given in [3]. Loosely speaking, (3.6) states that if there are no jump discontinuities occurring in the set [S.sub.m], then in regions of smoothness (away from the jumps) we have [L.sub.m]f([theta]) [right arrow] 0 with mth order accuracy. At the jump discontinuities (3.3) converges to the value of the jump.

As described in [3], the one-dimensional mth order edge detection method can be simplified to building the mth degree Newton's divided difference. Recall that for the local set defined in (3.2), the mth degree Newton divided difference for [S.sub.m] when m [greater than or equal to] 1 is defined by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3.7)

where f[[theta]] [equivalent to] f([theta]). We can compute (3.7) as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3.8)

for

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3.9)

As shown in [3], the edge detection method (3.3) can be expressed using Newton divided differences as

[L.sub.m]f([theta]) = [m!/[q.sub.m]([theta])] f[[S.sub.m]]. (3.10)

Furthermore, an explicit formula for the coefficients [c.sub.i]([theta]) can be derived as

[c.sub.i]([theta]) = m!/[[omega].sub.i]([S.sub.m]), i = 1, ... , m + 1. (3.11)

yielding the values [q.sub.m]([theta]) in (3.5). Details are found in [3].

Consider the following example displayed in Figure 3.1 sampled on 64 Gaussian points. Since the data are sparsely spaced, it is difficult to determine the steep gradients from the true jump discontinuities.

[FIGURE 3.1 OMITTED]

Example 3.1.

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3.12)

The results of the edge detection method are exhibited in Figure 3.2, where it is clear that [L.sub.m]f([theta]) [right arrow] 0 away from the discontinuities, and as predicted, the convergence increases with the order m.

However the application of (3.3) also encounters some problems in the approximation of jump functions. Specifically ms m increases, oscillations that occur in the neighborhood of a jump discontinuity can be misidentified as true edges, as is evident in Figure 3.2(d). On the other hand, for smaller m there is a risk of identifying a steep gradient as an edge, as seen in Figure 3.2(a). These problems can be exacerbated when the true function f([theta]) is not known and the reconstruction [f.sub.N]([theta]) must be used instead.

To avoid the possibility of misidentification of an edge, we adopt the minmod algorithm originally described in [17] and adapted in [3] for the local edge detection method. Jump discontinuities are distinguished from neighborhood oscillations by effectively incorporating information provided by the results from (3.3) using the various orders of m. This technique is referred to in [3] as the minmod local edge detection method.

Definition 3.1. For a given finite set M [subset] of N of positive integers, consider the set [L.sub.M] f = {[L.sub.m] f : R [right arrow] R | m [member of] M}. The mimnod function is defined by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (3.13)

Loosely speaking, the minmod function distinguishes artificial oscillations from existing jump discontinuities by recognizing that a true jump will not change sign for arbitrary m, while artificial oscillations will vary in sign depending on the order m chosen in (3.3). The application of the minmod function

Figure 3.3: The minmod algorithm (3.13) applied to Example 3.1 using M = 5. The original function is sampled on 64 Gaussian spaced grid points.

[FIGURE 3.2 OMITTED]

[FIGURE 3.3 OMITTED]

(3.13) is exhibited in Figure 3.3.

Even after applying (3.13), it is possible that some small oscillations may still persist in the jump function approximation. This is easily rectified by applying a critical threshold parameter that is scaled according to the particular image being analyzed as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3.14)

Here [L.sub.M]f defines the final jump approximation of f([theta]).

We now review the total procedure of the minmod local edge detection method. First, (3.3) is applied for various m. Then (3.13) and (3.14) are used to ensure that the jumps are correctly identified. The edges are then defined as each value [[theta].sub.i] such that [L.sub.M]f([[theta].sub.i]) [not equal to] 0 for i = 1, ... , [N.sub.[theta]]. The threshold value, [T.sub.crit], is problem dependent. However, [T.sub.crit] [approximately equal to] [DELTA][theta] is a reasonable choice assuming that the image is sufficiently resolved for the purpose of determining jump discontinuities. Specifically, we assume that only one jump discontinuity can exist between two neighboring points. If an error is made such that more edges are found than actually exist, then a high order reconstruction method will still work as the function is still smooth in each region of approximation. On the other hand, if an edge goes undetected, then the function will be approximated over a discontinuity resulting in the Gibbs phenomenon. In this case aposteriori error estimates should be used to determine if the edges must be recomputed.

To determine the edges of a piecewise smooth function from spherical harmonic data, we assume that we are given the grid point values f([[theta].sub.i], [[phi].sub.j]) for i = 1, ... , [N.sub.[theta]] and j = 1, ... , [N.sub.[phi]] over the domain [0, [pi]] x [0, 2[pi]]. One dimension is held fixed, e.g., [phi] = [[phi].sub.f], and then the edges are determined for each f([theta], [[phi].sub.f]) as the set

[[THETA].sup.*] = {[[theta].sup.*] : 0 [less than or equal to] [[theta].sup.*] [less than or equal to] [pi]}.

Each edge [[theta].sup.*] then constitutes an end point for a smooth interval of f([theta], [[phi].sub.f]). Specifically if

[[THETA].sup.*] = {[[theta].sup.+], [[theta].sup.*.sub.2], ..., [[theta].sup.*.sub.j]},

then

(0, [[theta].sup.*.sub.1]), ([[theta].sup.*.sub.1], [[theta].sup.*.sub.2]), ([[theta].sup.*.sub.2], [[theta].sup.*.sub.3]), ..., ([[theta].sup.*.sub.j], [pi])

form the smooth subintervals in which f([theta],[[phi].sub.f]) is to be reconstructed. The method is performed in exactly the same way in the longitudinal direction by holding [theta] = [[theta].sub.f] fixed. In this case, the edge detection algorithm is further simplified since the longitudinal data points [[phi].sub.j], j = 0, ..., [N.sub.[phi]] are uniformly distributed [3]. Hence, for any point ([theta], [phi]) in the domain [0, [pi]] x [0, 2[pi]], we can determine the smooth region which encloses it [[[theta].sub.a], [[theta].sub.b]] x [[[phi].sub.a], [[phi].sub.b]]. These will be the intervals used for the Gegenbauer reconstruction described in [section]4.

4 The Gegenbauer Reconstruction Method

In this section we review the basic concepts of the Gegenbauer reconstruction method [19, 20, 21]. Although we will ultimately be considering the reconstruction of a piecewise smooth function from spherical harmonics, it is useful to review the one-dimensional Gegenbauer reconstruction method for both trigonometric and general orthogonal polynomials expansions. Since the spherical harmonic partial stun expansion contains both Fourier and associated Legendre polynomials, the Gegenbauer reconstruction method must be amenable to both types of expansions.

As discussed in [section]1, the Gegenbauer reconstruction method can only be performed in regions of smoothness. Crossing over a discontinuity in a piecewise smooth function will lead to the Gibbs phenomenon. The endpoints of each smooth sub-region for both the latitudinal and longitudinal directions are first determined by the local edge detection method explained in [section]3.

4.1 Review of the Gegenbauer Polynomials

Before we describe the Gegenbauer reconstruction method, it is useful to briefly review the Gegenbauer polynomials and some related properties. More details can be found in [1, 5, 7, 18].

The Gegenbauer polynomials [C.sup.[lambda].sub.l]([xi]) on the interval [-1, 1] are orthogonal under weight function [(1 - [[xi].sup.2]).sup.[[lambda] - 1/2] for [lambda] [greater than or equal to] 0 such that

[[integral].sup.1.sub.1][(1 - [[xi].sup.2]).sup.[lambda]- 1/2][C.sup.[lambda].sub.k]([xi])[C.sup.[lambda].sub.l]([xi])d[xi] = [[delta].sub.l,k][h.sup.[lambda].sub.l],

where

[h.sup.[lambda].sub.l] = [[pi].sup.1/2][C.sup.[lambda].sub.l](1) [GAMMA]([lambda] + 1/2)/[GAMMA]([lambda])(l + [lambda]) and [C.sup.[lambda].sub.l](1) = [GAMMA](l + 2[lambda])/l![GAMMA](2[lambda]).

A function f([xi]) on [-1, 1] can be represented as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

As with other orthogonal polynomial expansions, the truncated approximation

[f.sub.m]([xi]) = [m.summation over (l=0)] [[??].sup.[lambda].sub.l][C.sup.[lambda].sub.l]([xi])

converges exponentially to f([xi]) as long as f([xi]) is smooth [7, 18]. Note that [lambda] - 1/2 and [lambda] = 0 yield the standard Legendre and Chebyshev polynomial expansions, respectively.

To reconstruct f(x) on a sub-interval x [member of] [a, b] we apply the simple linear transformation

x = [epsilon][xi] + [delta] (4.1)

where [epsilon] = b-a/2 and [delta] = b+a/2. In this case the Gegenbauer coefficients,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.2)

yield the approximation

[f.sub.m](x([xi])) = [m.summation over (l=0)][[??].sup.[lambda].sub.l,[epsilon]][C.sup.[lambda].sub.l]([xi]), (4.3)

which exponentially converges to f(x) in the interval [a, b] as m [right arrow] [infinity] when f(x) is smooth. In what follows, we review how the Gegenbauer polynomials can be applied to reconstruct piecewise analytic functions when the given data. prohibits the direct use of (4.2) and (4.3).

4.2 Development of the Gegenbauer Reconstruction Method

Assume that a function f([xi]) is piecewise smooth on an interval [-1,1] and analytic in the region a [less than or equal to] x [less than or equal to] b, where -1 [less than or equal to] a < b [less than or equal to] 1. We wish to reconstruct f(x) in the sub-interval [a, b] given either the Fourier coefficients,

[[??].sub.k] = 1/2 [[integral].sup.1.sub.-1] f([xi])[e.sup.-ik[pi][xi]] d[xi], k [member of] [- N/2, N/2] (4.4)

or orthogonal polynomial coefficients (e.g., Chebyshev or Legendre),

[a.sub.l] = 1/[gamma]l] [[integral].sup.1.sub.-1] f([xi])[P.sub.l]([xi])w([xi])d[xi], l [member of] [0, N], (4.5)

where [P.sub.l]([xi]) are orthogonal polynomials on [-1, 1] with weight function w([xi]) and normalization factor [[gamma].sub.l].

In the Fourier case, the approximation to f(x([xi])) using the transformation (4.1) is given by the partial sum

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.6)

Analogously,

[f.sub.N] (x([xi])) = [N.summation over (l=0)] [a.sub.l][P.sub.l](x([xi])) (4.7)

represents the standard orthogonal polynomial approximation. It is well known that if f([xi]) is smooth (and periodic in the Fourier expansion case) on [-1, 1], then the approximations (4.6) and (4.7) converge spectrally to f(x([xi])) for any interval [a, b]. However, if f is only piecewise smooth, the above approximations yield very poor results with Gibbs oscillations visible near the endpoints a and b as well as an overall reduced convergence throughout the interval [7, 18]. This is due to the fact that, even though the function may be smooth in the interval of reconstruction [a, b], the coefficients (4.4) and (4.5) are obtained by integrating a piecewise smooth function across the entire region [-1, 1], including the discontinuities. Furthermore, the Fourier approximation poses an additional problem since there is no reason to assume that f is periodic in [a, b].

The Gegenbauer reconstruction method, originally developed in [21], addresses this issue by observing that it is possible to reconstruct a function f(x) by reprojecting the approximations (4.6) and (4.7) onto the orthogonal Gegenbauer polynomials to construct (4.3) in a smooth sub-interval [a, b]. Clearly, the coefficients (4.2) are not known for this reconstruction; rather, (4.6) and (4.7) must be used to construct the approximate coefficients

[[??].sup.[lambda].sub.l,[epsilon]] = 1/[h.sup.[lambda].sub.l] [[integral].sup.1.sub.-1] [f.sub.N](x([xi]))[C.sup.[lambda].sub.l]([xi])[(1 - [[xi].sup.2]).sup.[lambda] - 1/2] d[xi] [approximately equal to] [[??].sup.[lambda].sub.l,[epsilon]], l [member of] [0, m]. (4.8)

Note that if (4.8) is derived from the Fourier sum (4.6), then there is an explicit formula [5] given by:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.9)

where [J.sub.l+[lambda]]([epsilon]k[pi]) is the Bessel function of the first kind [1]. If instead (4.8) is constructed from (4.7), we use the Chebyshev Gauss Labotto quadrature rule,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.10)

where the quadrature points are defined as ([[xi].sub.j] = -cos [[pi]j/N. For the coefficients (4.8) we have

q(x([[xi].sub.j])) = [f.sub.N]([epsilon][[xi].sub.j] + [delta])[(1 - [[xi].sup.2.sub.j]).sup.[lambda][C.sup.[lambda].sub.l]([[xi].sub.j])

and [??] [greater than or equal to] 2[lambda]+m+N. Another useful identity is the closed form for the Christoffel Darboux sum, which will be practical in numerically implementing (4.12) if the coefficients (4.8) are obtained from the orthogonal polynomial coefficients (4.5). That is:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.11)

where

[k.sub.m] = [2.sup.m][GAMMA]([lambda] + m)/m![GAMMA]([lambda]).

Once (4.8) is computed, the Gegenbauer reconstruction method

[g.sub.m](x([xi])) = [m.summation over (l=0)] [[??].sup.[lambda].sub.l,[epsilon]][C.sup.[lambda].sub.l]([xi]) (4.12)

approximates (4.3). As shown in [19, 20], (4.12) is an exponentially convergent approximation to f(x) in [a, b] provided A and m grow with N. (Note that this is why neither the Cheybshev or Legendre expansions in (4.7) can be used for the reprojection (4.12).) In this way, the Gibbs phenomenon can be completely removed up to the endpoints of each smooth sub-interval. Hence, the Gegenbauer method is much more powerful than standard filtering methods, since no blurring will occur at the "edges" of each region. In the investigations of [12, 15] the method has been shown to be effective for two-dimensional and spherical problems, as well as for various applications.

Remarks.

1. The method also applies when the given information are the pseudospectral coefficients rather than (4.4) or (4.5). This is particularly useful when the underlying image is given on grid-point data, as is the case for atmospheric spectral transform models.

2. Although theoretically m and [lambda] grow with N, in practice it has been shown that much smaller values will yield high resolution reconstructions. Smaller values of m and [lambda] are also necessary to avoid round off error and to reduce computational cost, [10, 11].

3. The Gegenbauer reconstruction method works well when the information is given as the partial sum expansion of general orthogonal polynomials (4.7), including associated Legendre polynomials, making it amenable to spherical harmonics, [12, 15].

We now wish to adapt the Gegenbauer reconstruction method to reconstruct piecewise smooth functions on spheres. This will be accomplished by reconstructing the grid-point values f([[theta].sub.i], [[phi].sub.j]), i = 1, ..., [N.sub.[theta]] and j = 1, ... , [N.sub.[phi]], over the whole domain [0, [pi]] x [0, 2[pi]]. Typically, the number of resolved grid-points, [N.sub.[theta]] and [N.sub.[phi]], are related to the number of given spherical harmonic coefficients (2.2). To demonstrate the Gegenbauer reconstruction method, we look a single point ([theta], [phi]) for which we want to reconstruct the function f. Let us assume that f([theta], [phi]) is analytic in the region [[theta].sub.a] [less than or equal to] [theta] [less than or equal to] [[theta].sub.b] and [[[phi].sub.a] [less than or equal to] [phi] [less than or equal to] [[phi].sub.b] which has been previously determined in [section]3.

To apply the Gegenbauer reconstruction method from spherical harmonics, we first make the linear transformations

[theta]([[xi].sub.[theta]]) = [[epsilon].sub.[theta]][[xi].sub.[theta]] + [[delta].sub.[theta]] [phi]([[xi].sub.[phi]]) = [[epsilon].sub.[phi]][[xi].sub.[phi]] + [[delta].sub.[phi]] (4.13)

where [[epsilon].sub.[theta]] = [[theta].sub.b]-[[theta].sub.a]/2, [[epsilon].sub.[phi]] = [[theta].sub.] - [[theta].sub.a]/2, [[delta].sub.[theta]] = [[theta].sub.b] + [[theta].sub.a]/2, and [[delta].sub.[phi]] = [[phi].su.b] + [[phi].sub.a]/2. The spherical harmonic partial sum in the region of reconstruction is

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.14)

As discussed in the introduction and visible in Figure 2.1, the spherical harmonic sum approximation (4.14) yields poor results, with Gibbs ringing particularly prevalent near the edges of the regions. However, noting that the spherical harmonic partial sum is a combination of Fourier and orthogonal polynomial expansions, we can apply the Gegenbauer reconstruction method in the same way as discussed above. Specifically, we use (4.14) to construct the coefficients

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.15)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

This leads to the spectrally convergent approximation [15]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.16)

4.3 Implementation of the Gegenbauer Reconstruction Method

In order for the Gegenbauer reconstruction method to be as effective as a postprocessor for geophysical and atmospherics simulations, it is important that the method be numerically efficient, robust, and easy to implement. These issues were first addressed in [11] and the resulting algorithm is briefly reviewed here.

To compute (4.15) efficiently, we first rewrite

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Applying (4.9) yields

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.17)

The Gauss Labotto quadrature formula is used to approximate

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.18)

where [c.sub.j] is defined in (4.10) for [xi][theta],j = -cos[pi]j/N, and [bar.N] [greater than or equal to] 2[[lambda].sub.1] + [m.sub.1] + N. Employing (4.17) and (4.18) enables us to rewrite (4.16) as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.19)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.20)

The direct application of the Christoffel Darboux formula (4.11) yields an explicit formula for the sum

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.21)

Hence, by defining

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.22)

we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.23)

There are several numerical advantages to computing (4.16) as (4.23). First, the Fourier expansion component of spherical harmonics (2.1) means that the explicit formula (4.17) can be applied instead of costly integration, which allows the use of FFT algorithms for the reconstruction in the longitudinal direction. Second, although we are forced to use a quadrature rule (4.10) to compute the integral in the latitudinal direction, the Christoffel Darboux formula (4.11) allows us to eliminate one of the sums in (4.20) which also reduces the cost of computation. Finally, by "splitting" the reconstruction (4.23) into latitudinal and longitudinal approximations, it is clear that the reconstruction can be performed by either fixing the latitudinal coordinate and computing

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.24)

or alternatively by fixing the longitudinal coordinate and computing

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.25)

Hence, the Gegenbauer reconstruction method can be tailored to a particular application. For example, if the global data field has more discontinuities in the longitudinal direction, it is reasonable to hold [theta] fixed and apply (4.24), and similarly for the latitudinal direction using (4.25). Clearly either of these approximations are more cost efficient than (4.23). We also note that it is possible to apply one of the reconstruction algorithms developed in [8, 9, 29] in the longitudinal direction. This study does not apply such a "splitting" technique to reconstruct a function, but this option will be explored in future investigations.

5 Numerical Examples

To demonstrate the effectiveness of the Gegenbauer reconstruction method for spherical harmonics, we recall Example 2.1. After determining the regions of smoothness by applying the local edge detection method in both the longitudinal and latitudinal directions, the Gegenbauer reconstruction method (4.16) is used to obtain a high-order reconstruction of the image without compromising the finer features of the function. Figure 5.1 shows the results. Note that the parameter choices [[lambda].sub.1] = [[lambda].sub.2] = 4, and [m.sub.1] = [m.sub.2] = 8 are not optimized, but rather chosen for ease of implementation while still ensuring convergence. Parameter optimization and the reduction of possible round off error for the Gegenbauer reconstruction method are discussed in [10, 14, 23]. As an example closer to real applications, we look at a picture of the surface mountains of the Earth. In this ease, although there are no actual discontinuities in the image, the resolution used is insufficient to properly resolve the steep gradients. Hence, as displayed in Figure 5.2(b), the Gibbs oscillations are visible in the model mountains field as a series of altitude oscillations over the oceans. They are especially prominent in the vicinity of sharp mountain ridges like the Andes. As discussed in [section]1, the Gibbs oscillations ill the mountain field is one of the causes of the instabilities arising in the atmospheric model coupled to the ocean model. Although filtering alleviates some of the ringing effects, it is clear from Figure 5.3(a) that the regions near the sharp mountain ridges are blurry, and some oscillations remain. The Gegenbauer reconstruction method (4.16), based on the spherical harmonic coefficients (2.2) for its reconstruction, is exhibited in Figure (5.3), where it is clear that the Gibbs oscillations are removed and accuracy is restored.

[FIGURE 5.1 OMITTED]

[FIGURE 5.2 OMITTED]

Remarks. There are several factors to be considered here.

[FIGURE 5.3 OMITTED]

1. The Gegenbauer reconstruction method should not be used in very small regions that contain too few grid points, as resolution requirements prohibit spectral reconstruction in this case. Local interpolation methods or splines could be used (with diminished accuracy in those regions) and will be discussed in future work. In these examples we used the filtered reconstruction values in the small regions.

2. As discussed in [section]4, it is also possible that one-dimensional reconstruction in either the latitudinal or longitudinal direction would yield sufficient accuracy. This would make the Gegenbauer reconstruction method less costly, particularly if the reconstruction was performed for longitudinal values since the FFT can be applied. Aposteriori error estimates can determine whether or not the full Gegenbauer reconstruction method should be applied, and will be discussed in future work.

3. No attempts at optimization were made, nor was a check of round off error enforced. This has been done in other applications [10, 14, 23], and could be useful here as well.

4. We have not attempted our algorithms for problems that have irregularities near the poles. This can cause added complications, since the reprojection method requires large enough smooth regions (with at least [pi] points per interval) in order to do the reconstruction. A similar issue was addressed in [2]. There the Gegenbauer reconstruction method was used to recover MRI brain images from Fourier data. In small smooth regions, where there were too few points to construct the Gegenbauer reconstruction, a

simple line segment was drawn across the interval. Perhaps something more sophisticated can be used on the sphere when the singularities occur near the pole.

6 Concluding Remarks

This study combines the edge detection algorithm [3], the minmod procedure [17], and the Gegenbauer reconstruction method, e.g., [19, 21, 12], to produce an automated high order reconstruction method for piecewise smooth functions using the spherical harmonics partial sum approximation (2.3). The numerical examples provided in [section] 3 and [section] 5 demonstrate the effectiveness of the method. By incorporating numerical techniques such as Newton divided differencing, FFT, Gaussian quadrature, and the Christoffel Darboux sum, tim method can be made numerically robust and cost efficient. Future work will include the application of the reconstruction method to real geophysical and climatology data, where it will be necessary to develop techniques that validate the reconstructed solutions. In addition, the method will be incorporated within the simulation of models.

References

[1] M. Abramovitz and I. Stegun, Handbook of Mathematical Functions, Dover Publications, New York, 1965.

[2] R. Archibald and A. Gelb, A Method to Reduce the Gibbs Ringing in MRI Scans While Keeping Tissue Boundary Integrity, IEEE Transactions on Medical Imaging, 21:4, 305-319, 2002.

[3] R. Archibald, A. Gelb, and J. Yoon, Polynomial Fitting for Edge Detection in irregularly Sampled Signals and Inmges, SIAM Journal on Numerical Analysis, 43, 259-279, 2005.

[4] R. Archibald, A. Gelb, and J. Yoon, Determining the Locations of the Discontinuities in the Derivatives of Functions, in press, 2006.

[5] H. Bateman, Higher Transcendental Functions, Vol. 2, McGraw-Hill, New York 1953.

[6] J.P. Boyd, Chebyshev and Fourier Spectral Methods, Springer-Verlag, 1989.

[7] C. Canuto, M.Y. Hussaini, A. Quarteroni, and T.A. Zang, Spectral Methods in Fluid Dynamics, Springer 1987.

[8] T.A. Driscoll and B. Fornberg, A Pade-based Algorithm for Overcoming the Gibbs Phenomenon, emphNumerical Algortihms, 26, 77-92, 2001.

[9] K. S. Eckhoff, On a High Order Numerical Method for Functions with Singularities, Mathematics of Computation, 67, 1063-1087, 1995.

[10] A. Gelb, Parameter Optimization and Reduction of Round Off Error for the Gegenbauer Reconstruction Method, Journal of Scientific Computing, 20:3, 433-459, 2004.

[11] A. Gelb, A Hybrid Approach to Spectral Reconstruction of Piecewise Smooth Functions, Journal of Scientific Computing, 15, 293-322, 2001.

[12] A. Gelb, The Resolution of Gibbs Phenomenon for Spherical Harmonics, Mathematics of Computation 66:218, 699-717, 1997.

[13] A. Gelb and J. Gleeson, Spectral Viscosity for Shallow Water Equations in Spherical Geometry, Monthly Weather Review, 129:9, 2346-2360, 2001.

[14] A. Gelb and Z. Jackiewicz, Determining Analyticity for Parameter Optimization of the Gegenbauer Reconstruction Method, SIAM J. Sci. Comput 27:3, 1014-1032, 2006.

[15] A. Gelb and A. Navarra, Recovering Grid-Point Values without Gibbs Oscillations in Two Dimensional Domains on the Spheres, Center for Research on Parallel Computation, California Institute of Technology, 1997.

[16] A. Gelb, and E. Tadmor, Detection of Edges in Spectral Data II. Nonlinear Enhancement, SIAM J. Namer. Anal. 38, 1389-1408, 2000.

[17] A. Gelb, and E. Tadmor, Adaptive Edge Detectors for Piecewise Smooth Data Based on the MinMod Limiter, Journal of Scientific Computing, 28-23, 279-306, 2006.

[18] D. Gottlieb and S. Orszag, Numerical Analysis of Spectral Methods: Theory and Applications. SIAM, Philadelphia, 1977.

[19] D. Gottlieb and C.-W. Shu, A General Theory for the Resolution of the Gibbs Phenomenon, w, Tricomi's Ideas and Contemporary Applied Mathematics Convention, Roma Accademia Nazionale Dei Lincei, Italy, 1998.

[20] D. Gottlieb and C.-W. Shu, On the Gibbs Phenomenon and its Resolution, SIAM Review, 1997.

[21] D. Gottlieb, C.-W. Shu, A. Solomonoff, and H. Vandeven, On the Gibbs Phenomenon I: Recovering Exponential Accuracy from the Fourier Partial Sum of a Non-Periodic Analytic Function, J. Comput. Appl. Math., 43, 81-92, 1992.

[22] G.J. Haltiner and R.T. Williams, Numerical Predictions and Dynamics Meteorology, Wiley, 1980.

[23] Z. Jackiewicz, Determination of optimal parameters for the Chebyshev-Gegenbauer reconstruction method, SIAM J. Sci. Comput., 25:4, 11871198, 2004.

[24] C. Lindberg and A. Broccoli, Representation of topography in spectral climate models and its effect on simulated precipitation, J. Climate, 9, 26412659, 1996.

[25] Y. Maday, S.M. Ould Kaber, and E. Tadmor, Legendre Pseudospectral Viscosity Method for Nonlinear Conservation Laws, SIAM J. Num. Anal., 30, 321-342, 1993.

[26] F.B. Mitchell, The Greenhouse Effect and Climate Change, Rev. of Geophys., 27, 115-139, 1989.

[27] A. Navarra, W.F. Stern, and K. Miyakoda, Reduction of the Gibbs Oscillations in Spectral Model Simulations, J. Climate, 7, 1169-1183, 1994.

[28] S.A. Orszag, Fourier Series on Spheres, Mon. Weather Rev. 102:56, 1974.

[29] B.D. Shizgal and J.H. Jung, Towards the Resolution of the Gibbs Phenomena, Journal of Computational and Applied Mathematics, 161, 41-65, 2003.

[30] C.W. Shu and P. Wong, A note on the accuracy of spectral method applied to nonlinear conservation laws, Journal of Scientific Computing, 10, 357369, (1995).

[31] P.N. Swarztrauber, On the Spectral Approximation of Discrete Scalar and Vector Functions on a Sphere, Siam J. Numer. Anal., 16, 934-949, 1979.

[32] E. Tadmor, Super Viscosity and Spectral Approximations of Nonlinear Conservation Laws, Numerical Methods for Fluid Dynamics IV, M. J. Baines and K. W. Morton eds., Clarendon Press, Oxford, 1993.

Chris Blakely

The University of Maryland, CSCAMM

4146 CSIC Building # 406, Paint Branch Drive, College Park, MD 20742-3289

cblakely@cscamm.umd.edu

Anne Gelb

Department of Mathematics and Statistics, Arizona State University

Tempe, Arizona 85287

ag@math.la.asu.edu

Antonio Navarra

National Institute of Geophysics and Volcanology (INGV)

Via Donato Creti 12. [-40128 Bologna, Italy

navarra@bo.ingv.it

Printer friendly Cite/link Email Feedback | |

Author: | Blakely, Chris; Gelb, Anne; Navarra, Antonio |
---|---|

Publication: | Sampling Theory in Signal and Image Processing |

Article Type: | Report |

Geographic Code: | 1USA |

Date: | Sep 1, 2007 |

Words: | 7512 |

Previous Article: | Modeling of power electronics circuits using wavelet theory. |

Next Article: | Tight frame completions with prescribed norms. |

Topics: |