Georges-Pierre Bonneau

Optimal Triangular Haar Bases for Spherical Data

Keywords: visualization, multiresolution, wavelets, orthogonality

Multiresolution analysis based on FWT (Fast Wavelet Transform) is now widely used in Scientific Visualization. Spherical biorthogonal wavelets for spherical triangular grids where introduced in [5]. In order to improve on the orthogonality of the wavelets, the concept of nearly orthogonality, and two new piecewise-constant (Haar) bases were introduced in [4]. In our paper, we extend the results of [4]. First we give two one-parameter families of triangular Haar wavelet bases that are nearly orthogonal in the sense of [4].

Then we introduce a measure of orthogonality. This measure vanishes for orthogonal bases. Eventually, we show that we can find an optimal parameter of our wavelet families, for which the measure of orthogonality is minimized. Several numerical and visual examples for a spherical topographic data set illustrates our results.

INTRODUCTION

Representing large data sets by decomposition in a wavelet basis has many applications in scientific visualization. This includes LOD (level of detail) visualization, local area zooming, progressive transmission over networks, compression, and others ([6]). Since wavelets were originally discovered in signal analysis, where the data comes from regular samples, their application in scientific visualization began with uniformly gridded data sets. An important application domain for uniformly gridded data visualization is volume visualization ( [START_REF] Gross | Two methods for Wavelet-Based volume rendering[END_REF]).

In the last few years, a major challenge has been to extend wavelet decompositions to work on more general data sets. The introduction of triangular wavelets in [START_REF] Lounsbery | Multiresolution analysis for surfaces of arbitrary topological type[END_REF] for geometric meshes, and triangular spherical wavelets in [START_REF] Schröder | Spherical wavelets: Efficiently representing functions on the sphere[END_REF] for scalar spherical data sets, were two crucial steps in that challenge. More recently, generalized Haar wavelets for irregular triangulations have been introduced in [START_REF] Bonneau | Multiresolution analysis on irregular surface meshes[END_REF].

Triangular wavelets are defined over nested triangular grids. These grids are the result of a subdivision process starting on a base triangular mesh, and dividing recursively each triangle into four subtriangles by inserting new vertices at the midpoint of each edge. In LMC-CNRS, BP.53, F-38041 Grenoble Cedex 9, France.

the case of spherical data, each edge is part of a great circle, the midpoint is the geodesic bisector, and the base mesh may be either a spherical tetrahedron, octahedron or icosahedron, see Fig. 1.

2 subdivision steps starting on a spherical tetrahedron 2 subdivision steps starting on a spherical octahedron 2 subdivision steps starting on a spherical icosahedron Figure 1: Subdivision of spherical triangular meshes This paper deals with piecewise constant triangular wavelets. While this may be seen as a restriction for smooth data sets, it is sufficient for large and non-smooth data. In particular, many earth data sets, which are of special interest in scientific visualization, are usually non-smooth, and thus can be tackled efficiently with piecewise constant wavelets. This has been verified in [START_REF] Schröder | Spherical wavelets: Efficiently representing functions on the sphere[END_REF] where Haar bases are shown to work as good as smoother bases for the compression of the topographic data set that we will use in our illustrations.

Wavelet bases may be orthogonal, semi-orthogonal or simply biorthogonal. As we shall see later in more detail, over these three types, orthogonal wavelet is the only one which gives direct access to best approximation by discarding wavelet coefficients. Unfortunately, there is no orthogonal triangular Haar bases on the sphere. To overcome this lack, the concept of nearly orthogonality was introduced in [START_REF] Nielson | Haarwavelets over triangular domains with applications to multiresolution models for flow over a sphere[END_REF]. A spherical wavelet basis is said to be nearly orthogonal if it becomes orthogonal when the subdivision depth increases (i.e. when the spherical triangles become planar). Two wavelet bases that are nearly orthogonal were given in [START_REF] Nielson | Haarwavelets over triangular domains with applications to multiresolution models for flow over a sphere[END_REF].

A u t h o r p r e p r i n t

In this paper, we extend the results of [START_REF] Nielson | Haarwavelets over triangular domains with applications to multiresolution models for flow over a sphere[END_REF] in the following directions. First, we give two one-parameter families of spherical Haar wavelet bases that are nearly orthogonal and simpler than in [START_REF] Nielson | Haarwavelets over triangular domains with applications to multiresolution models for flow over a sphere[END_REF]. Then we introduce a measure of orthogonality. This measure vanishes for orthogonal bases, and converges to zero for nearly orthogonal bases, as the subdivision depth increases. Eventually, we show that there exist an optimal parameter value for our new wavelet families, for which the measure of orthogonality is minimal. This paper is structured as follows. In section 2, we give the basics on piecewise constant triangular Haar wavelets. Section 3 is a review of the concept of orthogonality, semi-orthogonality, and nearly orthogonality in the context of triangular Haar wavelets. In section 4, the two new nearly orthogonal one-parameter families of wavelets are introduced. Section 5 is dedicated to the measure of orthogonality and to the optimization of this measure over the new families of wavelets. Eventually, section 6 illustrates our results on an earth topography data set.

Piecewise constant triangular wavelets

Piecewise constant triangular wavelets are defined by their values on the four subtriangles dividing their support. For each triangle in the nested grid, there exist three triangular Haar wavelets whose support matches exactly with that triangle. The global decomposition/reconstruction process of a data set in this wavelet basis consists in the recursive application of a local decomposition/reconstruction process.

Local decomposition/reconstruction

Let be a triangle at the subdivision depth ( for the base mesh). Let , , , be the four sub-triangles of . Let , and denote the three wavelet functions whose support matches , and denote the value of on the sub-triangle . Eventually, let denote the constant function whose value is 1 on . Fig. 2 illustrates these notations, and section 3 includes a practical numerical example of wavelet functions , and

. The fundamental property of triangular Haar wavelets is that every piecewise constant function on the four sub-triangles ( , ) can be written as a linear combination of the constant function and the three wavelets:

x x x x 1 2 3 0 k+1 k+1 k+1 k+1
The relation that gives the values of the left-hand side function from the wavelet coefficients is referred to as local reconstruction:

(1)

Let denote the so-called reconstruction matrix in [START_REF] Bonneau | Multiresolution analysis on irregular surface meshes[END_REF]. The inverse process of ( 1) is referred to as local decomposition:

(2)

Global decomposition / reconstruction

Given the values of an input data set at the maximal depth level , the decomposition in the wavelet basis can be described by the following pseudo-code:

Orthogonality, semi-orthogonality and nearly orthogonality

We will review the concepts of orthogonality, semi-orthogonality and nearly orthogonality in the context of spherical triangular Haar wavelets, although the first two of them are much more general.

Orthogonality / semi-orthogonality

Here we use the notations of section 2. A function at level is decomposed in the wavelet basis:

x x x x 1 2 3 0 k+1 k+1 k+1 k+1
The wavelet basis is defined as semi-orthogonal, if the inner product of the wavelets with the constant function vanishes:

The wavelet basis is defined as orthogonal, if in addition to the semiorthogonality condition, any two distinct wavelet functions have a vanishing inner product:

A u t h o r p r e p r i n t

We can illustrate the weakness of semi-orthogonal bases on a simple example. Assume the triangles are planar, and the areas of the four subtriangles equal identically one. In this case, the semi-orthogonal basis introduced in [START_REF] Schröder | Spherical wavelets: Efficiently representing functions on the sphere[END_REF] is the following:

1 1 1 -3 1 1 1 -3 1 1 1 -3 1
We wish to represent the following function using this basis:

1 3 0 0
The exact decomposition of in the wavelet basis is:

The approximation of after discarding the lowest wavelet coefficient (compression gain: 25%) is:

5/2 5/2 -1/2 -1/2
The L error for this approximation is .

The best approximation of with the same three basis functions is:

2 2 0 0
In this case, the L error is reduced to . Therefore, approximation by discarding the lowest wavelet coefficient, results in an L error that is larger than the best L error with the same wavelet functions. This illustrates one advantage of orthogonal bases over semiorthogonal bases: the best error obtained by keeping a certain subset of the wavelet functions is always obtained for the corresponding wavelet coefficients.

Nearly orthogonality

The concept of nearly orthogonality for spherical wavelets was introduced in [START_REF] Nielson | Haarwavelets over triangular domains with applications to multiresolution models for flow over a sphere[END_REF]. As the subdivision depth increases, the spherical triangles become planar, and the areas of the four sub-triangles become uniform. On the other hand, the three wavelet functions depend on the areas of these sub-triangles. A wavelet basis is defined as nearly orthogonal, if it is orthogonal in the limit case, i.e. when the four sub-triangle areas are equal.

Families of nearly orthogonal spherical Haar wavelets

In this section, we introduce new spherical Haar wavelets, that are nearly orthogonal, and that depend on one free parameter. In addition to the previous notations given in section 2, we denote , , , the areas of the four sub-triangles , , , .

Previous spherical triangular Haar bases

In order to simplify the notations, we will always give the unnormalized reconstruction matrices: to find the effective, L -normalized reconstruction matrices, each column has to be divided by the square rooted sum of the squared coefficients weighted by the areas; and each line has to be multiplied by . In other words, if is the unnormalized reconstruction matrix, the effective reconstruction matrix is given by:

In [START_REF] Schröder | Spherical wavelets: Efficiently representing functions on the sphere[END_REF], the following unnormalized reconstruction matrix was introduced:

(3)

This corresponds to a semi-orthogonal basis, since with this matrix.

This basis is not nearly orthogonal, since in the limit case of uniform areas, i.e. when , we have the following inner products:

In [START_REF] Nielson | Haarwavelets over triangular domains with applications to multiresolution models for flow over a sphere[END_REF], the following unnormalized reconstruction matrices were proposed:

(4) and [START_REF] Schröder | Spherical wavelets: Efficiently representing functions on the sphere[END_REF] where

. Both of these matrices were proved to be semi-orthogonal and also nearly orthogonal in [START_REF] Nielson | Haarwavelets over triangular domains with applications to multiresolution models for flow over a sphere[END_REF].

New families of nearly orthogonal Haar bases

A quick look at previously introduced reconstruction matrices (3), ( 4), [START_REF] Schröder | Spherical wavelets: Efficiently representing functions on the sphere[END_REF] shows the 4 3 sub-matrices containing the wavelet function values are polynomials in the sub-triangle areas . The degree is 2 for the matrices (4) and ( 5) and 1 for the matrix (3).

Since we wanted to find both simple (i.e. with low degree polynomial coefficients), and nearly orthogonal matrices, we have started with a 4 3 sub-matrix containing arbitrary polynomials of degree 1 (linear functions) in . Without any condition on the coefficients of these polynomials, the number of free parameters would be 4 3 4 = 48.

A first condition that has to be fulfilled is symmetry. This property states that, when permuting the indices (1,2,3) in the right and left-hand side of the local decomposition/reconstruction relations (1) and ( 2), the equality should be preserved. Imposing symmetry reduces the number of free parameters from 48 to 10.

The second condition is semi-orthogonality. This condition reduces further the number of free parameters to three, and we get the following unnormalized, semi-orthogonal families of reconstruction matrices:

(6)

Here are the free parameters. As a special case of (6), the semi-orthogonal reconstruction matrix (3) of [START_REF] Schröder | Spherical wavelets: Efficiently representing functions on the sphere[END_REF] is obtained for , and .

The last condition to impose is nearly-orthogonality. In the limit case of uniform areas, i.e. , the inner products of the wavelet functions yield:

These inner products vanish if and only if or Thus, we have two new families of symmetric, semi-orthogonal and nearly orthogonal triangular wavelets:

(7) (8)
It should be kept in mind that (7) and (8) must be normalized, and this implies that there is really only one free parameter in (7) and (8). In fact, we can assume without restriction that in the rest of this paper (more precisely, if

, then replacing by one, and by gives the same matrix after normalization). Thus ( 7) and (8) define two one parameter families of wavelet bases that are symmetric, semi-orthogonal, and orthogonal if .

Orthogonality measure and optimal Haar bases

In section 4, three previous wavelet bases were recalled, and two new families of wavelet bases were introduced. We now introduce a measure to distinguish all these wavelet bases and we show that, for the wavelet families ( 7) and ( 8), optimal values with respect to this measure can be calculated.

Since our main concern is on the orthogonality of the bases, the following measure was chosen:

(9)

Here is the finest subdivision level, is the number of triangles at that level ( equals to the number of base faces times 4 ), traverses all levels, and traverses all triangles at the level .

A straightforward property of the orthogonality measure ( 9) is that it vanishes if and only if all inner products also vanish, i.e. if the wavelet basis is orthogonal.

Another property of this measure is illustrated in Fig. 3. This figure shows the behaviour of the orthogonality measure (9) for the semiorthogonal basis (3) introduced in [START_REF] Schröder | Spherical wavelets: Efficiently representing functions on the sphere[END_REF], and for the new nearly orthogonal wavelets ( 7) with parameters and . The measure (9) has been computed at increasing subdivision depths , and for three base spherical triangulations (tetrahedron, octahedron, icosahedron). The behaviour of the measure (9) is the same for the three nested triangulations: for the semi-orthogonal basis (3) it increases up to a limit, and for the nearly orthogonal basis (7) with and it decreases and converges to zero. This is due to the fact that the inner products in (9) converge to a non-zero value for the semiorthogonal basis (3) and to zero for all nearly orthogonal bases.

Thus we see that the orthogonality measure (9) allows to distinguish the semi-orthogonal and the nearly orthogonal bases. It can also be observed that the measure decreases more quickly for the nested grid based on an icosahedron. This is because the triangles become more quickly uniform for this grid. . The measure (9) is computed at increasing subdivision depths and for three different nested spherical triangulations.

In the rest of this section, we use the orthogonality measure (9) to find an optimal basis between the two new families of bases (7) and (8).

In order to do this, we have fixed a certain subdivision depth , and compute the orthogonality measure (9) for different values of 1: These tables give the optimal parameter values for the two new families of wavelet bases ( 7) and ( 8), and the corresponding optimal orthogonality measure. The optimal parameter values have been computed for three different nested triangulations, at the subdivision depth 5, and up to a precision .

the free parameters in the bases ( 7) and (8). Fig. 4a shows a plot of the results for the family (7), with parameters and increasing from -50 to 50, for the three different nested triangulations, and at the subdivision depth 3. This plot suggests that there exist optimal parameter values for which the orthogonality measure (9) is minimal. Fig. 4b shows a zoom for the parameter values , increasing from -2 to -1. On this plot, the optimal parameter values for the family (7) are now clearly visible. The same behaviour is shown on figures 4c and 4d for the wavelet family (8).

Figure 4 shows also that the family (8) enables to reach a significantly smaller orthogonality measure. In table 1, the computed optimal values at the depth 5, up to a precision , together with the corresponding orthogonality measure (9) are given.

Results

The new nearly orthogonal bases (7) and ( 8) with optimal values given in table 1 have been tested and compared against the previous semi-orthogonal basis [START_REF] Lounsbery | Multiresolution analysis for surfaces of arbitrary topological type[END_REF]. The tests were performed on three nested spherical triangulations (base mesh tetrahedron, octhedron and icosahedron), subdivided to the level 7. The original topographic data set, available from the National Oceanographic and Atmospheric Administration 1 is a 4320 2160 uniform grid of short values. This data set was mapped on the finest triangulation in order to find the input data values of the wavelet decomposition algorithm. The mapping was done by a simple trapezoid quadrature formula that computes the integral of the original data on the finest triangles. The data set was then decomposed and reconstructed with an increasing number of wavelet coefficients. For each test, the difference between the true L error and the square rooted sum of the squared wavelet coefficients was computed. This measures the orthogonality of the wavelet basis with respect to the input data set in the following sense: it vanishes for orthogonal bases, and increases when the inner products of the wavelets modify the true L error. 2: These tables give the relative difference between the true error and the square rooted sum of the squared wavelet coefficients used for the reconstruction of the topographic data set. This difference vanishes for orthogonal bases, and should be smaller as possible. The tables show, in particular for the spherical icosahedron base mesh, that the difference is infinitesimal for our optimal wavelet bases (7) and (8).

wavelet bases ( 7) and (8) with the optimal parameters given in table 1 greatly reduces the L error due to the non-orthogonality of the basis functions.

The color plate shows some partial reconstructions obtained with the wavelet basis (8) and optimal parameter values, for the spherical icosahedron nested triangulation. The base mesh can be seen in red on each figure.

A u t h o r p r e p r i n t

Figure 2 :

 2 Figure 2: basis functons at level

  -orthogonal basis octahedron, nearly-orthogonal basis icosahedron, nearly-orthogonal basis tetrahedron, semi-orthogonal basis octahedron, semi-orthogonal basis icosahedron, semi-orthogonal basis

Figure 3 :

 3 Figure 3: Orthogonality measure (9) for the semi-orthogonal basis (3) introduced in [5] and for the new nearly-orthogonal wavelet basis (7) with parameters . The measure (9) is computed at increasing subdivision depths and for three different nested spherical triangulations.

[ 6 ]Figure 4 :

 64 Figure 4: Orthogonality measure for the new wavelet bases (7) (figs (a) and (b)) and (8) (figs (c) and (d)), computed at the subdivision depth 3.The parameter is fixed to , and the parameter is chosen on the X-axis. These plots show the existence of optimal parameter values for which the orthogonality measure is minimal.

Table 2 gives

 2

	# wav. coeff.	comp. gain	basis (3) basis (7) basis (8)
	500	99.2%	1.11 %	0.32 %	0.036 %
	2000	96.9 %	1.76 %	0.30 %	0.034 %
	6000	90.8 %	1.47 %	0.30 %	0.033 %
	10000 84.7 %	1.04 %	0.30 %	0.033 %
		Base mesh: spherical tetrahedron
	# wav. coeff.	comp. gain	basis (3) basis (7)	basis (8)
	1000	99.2 %	1.27 %	0.05 %	0.0056 %
	5000	96.2 %	1.36 %	0.05 %	0.0052 %
	10000 92.4 %	0.78 %	0.05 %	0.0052 %
	20000 84.7 %	0.58 %	0.05 %	0.0051 %
		Base mesh: spherical octahedron
	# wav. coeff.	comp. gain	basis (3)	basis (7)	basis (8)
	5000	98.5 %	2.44 %	0.0039 % 0.00058 %
	10000 96.9 %	2.30 %	0.0040 % 0.00054 %
	20000 94.0 %	1.90 %	0.0040 % 0.00052 %
	40000 87.8 %	1.59 %	0.0040 % 0.00051 %
		Base mesh: spherical icosahedron
	Table

the numerical results. It clearly shows that our new 1 at ftp://ftp.ngdc.noaa.gov/Solid Earth/Topography/tbase 5min/


G P Bonneau

Multiresolution Analysis on Irregular Surface Meshes

Keywords: wavelets, non-regular triangulations, compression, visualization

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

This paper is dedicated to the multiresolution analysis of piecewise constant data sets defined on irregular triangular surface meshes. It aims to give an efficient representation of such data sets, based on decomposition and reconstruction algorithms with respect to a specific hierarchical structure. Hierarchical structures play a crucial role since the beginning of computer graphics. More recently multiresolution analyses based on wavelet theory have started a new research domain closely related to previous works on hierarchical methods for computer graphics ( [START_REF] Sollnitz | Wavelets for Computer Graphics: Theory and Applications[END_REF]). Very briefly, the idea of wavelet techniques is to code a data set as a crude approximation, followed by a sequence of so-called detail coefficients that measure the error between successive finer approximations (see a survey on wavelet methods in [START_REF] Jawerth | An overview of wavelet based multiresolution analyses[END_REF]).

LMC-CNRS, B.P. 53, F-38041 Grenoble Cedex 9, France. http://www-lmc.imag.fr/˜bonneau/ Wavelet methods assume that the mesh on which the data is defined can be reached by recursive subdivision of a basic mesh. Thus every wavelet-based scheme is associated to hierarchies that have a tree structure (where every parent node is subdivided into a set of child nodes): wavelet volume visualization ( [START_REF] Gross | Two methods for Wavelet-Based volume rendering[END_REF]) is related to octree structures, wavelet radiosity ( [START_REF] Gortler | Wavelet radiosity[END_REF]) and wavelet over triangulated domains ([10, 12, 11]) are based on quadtree structures... On the other hand, irregular triangular meshes cannot be reached by subdivision rules, therefore hierarchical structures that have been developed to handle them are more complicated than trees. These include for example hierarchical Delaunay triangulations ( [START_REF] Kirkpatrick | Optimal search in planar subdivisions[END_REF][START_REF] De | A pyramidal data structure for trianglebased surface description[END_REF]), or progressive meshes ( [START_REF] Hoppe | Progressive meshes[END_REF][START_REF] Hoppe | View-dependent refinement of progressive meshes[END_REF]). This paper fills a gap between wavelet methods (on subdivision hierarchies) and hierarchical structures on irregular triangular meshes, for a special type of data set. In section 2 a wavelet-like transform is introduced, that takes as input a piecewise constant function on a triangulated domain, and outputs an approximation on a simpler triangulation of the same domain, together with a set of detail coefficients encoding the error with the input function. This decomposition reduces to a traditional Haar-wavelet transform, if the input triangulation subdivides the output triangulation. In section 3 this decomposition is used locally together with a hierarchical Delaunay triangulation, in order to define a global decomposition algorithm for piecewise constant data sets on irregular surface meshes. Two types of progressive exact reconstructions are illustrated on examples. Finally some improvements of the decomposition/reconstruction algorithms are described in section 4, and illustrated on "real size" data sets (1.3M and 920K faces).

The triangular wavelet-like transform

In this section we first explain why wavelet theory can't be applied for the multiresolution analysis of data sets on irregular surface triangulations. Then we explain how to extent wavelet theory in order to enable the analysis of piecewise constant data sets on irregular surface triangulations.

Need for an extension of wavelet theory

There exist numerous algorithms for reducing the number of triangles in a mesh -independently of the data that is defined on this mesh. Delaunay-removal can be applied to planar or convex triangular meshes, edge-collapse to general triangular meshes. If the mesh comes from the recursive 4-split of some triangles in a base mesh, then the obvious way to simplify it, is to replace each group of four subtriangles by their parent triangle. The common setting of these decimation algorithms is that a set of triangles is replaced by a set of triangles covering the same domain, with , as shown in fig. 1. If there exist data defined on the finest mesh, three main questions arise:

How should we compute the data on the new simplified meshes, in order to achieve good/best approximation of the original data ? A u t h o r p r e p r i n t Where should the triangle decimation be applied, in order to reach high compression ratio for a given error ?

Is it possible to recover exactly the original data set, from any approximation ?

Since recent works ( [START_REF] Lounsbery | Multiresolution analysis for surfaces of arbitrary topological type[END_REF][START_REF] Schröder | Spherical wavelets: Efficiently representing functions on the sphere[END_REF]), it is well known that wavelet theory can answer these three questions if the triangular mesh is built by recursive 4-split. The basic fact that enables to apply wavelet theory for these special meshes is that finer triangulations subdivide coarser ones. This is obviously not the case with Delaunay-removal or edgecollapse. Since we want to answer the three preceding questions for data sets defined on irregular surface meshes, we have to cope with the non subdivision connectivity of the meshes generated through the decimation algorithm. This is possible only through an extension of wavelet theory for which the nestedness of the sequence of approximating spaces is not required anymore. This extension was introduced in a previous paper [START_REF] Bonneau | BLaC-Wavelets: A multiresolution analysis with non-nested spaces[END_REF] in a different context. The reader should refer to [START_REF] Bonneau | BLaC-Wavelets: A multiresolution analysis with non-nested spaces[END_REF] for a context-independent presentation of this extension.

In the following, we apply this extension for the analysis of piecewise constant data sets on irregular surface meshes.

The local decomposition and reconstruction algorithm

First some terminology: the symbol denotes a triangle on the domain surface, denotes a data value (a scalar) on a triangle, denotes a detail coefficient. The superscript (fine) denotes quantities before the local decimation algorithm, and (coarse), after the decimation. Bold letters denote vectors. The pair denotes the piecewise constant function equal to on the triangles . If is a matrix, then denotes the -th column vector of . We focus here on the following setting: we are given a piecewise constant function , on triangles, and a set of triangles covering the same domain on the surface, with .

Analysis and synthesis formulas for irregular triangulations.

The essence of our local decomposition and reconstruction algorithm is the same as in wavelet theory: a function living in a fine space (in our case the piecewise constant function on the finer triangulation) is decomposed in a coarser approximating function (piecewise constant on the coarser triangulation) and error functions (piecewise constant on the finer triangulation). These error functions have two properties: they can be used to recover exactly the original data, and their norm is a measure of the error between the input function and the approximation.

Intuitively, we have to define one smoothing operator, that maps the input function on its approximation, and one error operator, that captures the difference between the input function and its approximation. These two operators are defined by two rectangular matrices and of size and respectively:

(1)

The smoothing operator (1) computes the coarser coefficients from the finer coefficients , and the error operator (2) computes the detail coefficients . The actual computation of the so-called analysis matrices and is detailed in the next section.

In order to keep a constant memory size for the data values, the original coefficients are cleared from memory after the decimation, and replaced by the coarse and detail coefficients and . Of course, the sum of the sizes of and equals the size of . Since is cleared from memory, the decomposition formulas (1) and ( 2) have to be invertible, in order to be able to recover the original data values. This is the purpose of the reconstruction formula:

(3)

The so-called synthesis matrices and are of sizes and respectively. Intuitively, the operator is the inverse of the smoothing operator : acts as a subdivision operator, although subdivision is not possible if the triangular domains are non-nested. The operator adds the details to the oversampled data , in order to recover the original data . The matrices and can be computed from and by: To be more precise about the properties of our decomposition, let us rewrite the reconstruction formula (3) with a functional point-ofview instead of a coefficient point-of-view: [START_REF] Gortler | Wavelet radiosity[END_REF] The next remarks can be made on (4), depending on whether the fine triangulation is a subdivision of the coarser one or not:

In any case, ( 4) is an orthogonal decomposition, and is always the best approximation of both and .

If

is a subdivision of , then , and (4) can be interpreted as a Haar-wavelet decomposition. In this case, the functions are the orthonormal Haar wavelet functions.

On the other hand, if is not a subdivision of , then the function can be interpreted as an intermediate approximation between the finer one and the coarser one . This intermediate function acts as a link between two functional spaces that are not nested. And the detail coefficients measure the error between this link and the finer approximation.

Before going further with the computation of the analysis and synthesis matrices, we insist on two important features of our local decomposition and reconstruction:

It can be used for every two sets of triangles , covering the same surface domain.

It leads to a Haar wavelet decomposition if

is a subdivision of .

Computation of the analysis and synthesis matri-

ces.

Computation of .

The analysis matrix in (1) has to be chosen so that is the best approximation of . This leads to the following usual choice in wavelet theory: must be the matrix of the scalar products between the dual coarser scaling functions and the primal finer scaling functions. (see for ex. chap. 7.4.1. in [START_REF] Sollnitz | Wavelets for Computer Graphics: Theory and Applications[END_REF] for the definition of primal and dual basis functions). In our case, the dual coarser scaling functions equal area on the triangle and vanish elsewhere. This implies the following analysis matrix : area area contains the relative areas of the intersections between fine and coarse triangles (see fig. 4). Therefore the value on a coarse triangle is the mean of the values on the fine triangles it intersects, weighted by the area of the intersection with each of such triangles. In the planar case, one careful way to compute is to use barycentric coordinates with respect to the vertices of , before computing the intersections. That way, all coordinates of the intersections lie between and , and no division has to be done to compute the relative area.

Computation of .

To ensure the orthogonality in the decomposition (3), the functions must be orthogonal to , for all pairs , with respect to the functional scalar product. In other words, the column vectors of must be orthogonal to the column vectors of , for the following scalar product: area It can be easily shown that this is equivalent to the lines of being orthogonal to the lines of for the following scalar product: area [START_REF] Gross | Two methods for Wavelet-Based volume rendering[END_REF] This leads to an underdetermined homogeneous linear system of equations. To solve it uniquely, we fix a submatrix of to the identity, and we inverse the resulting linear system. This yields a first choice for the lines of . In a second step, we perform a Gram-Schmidt orthogonalization process (see [START_REF] Strang | Introduction to Applied Mathematics[END_REF], Chap. 5.2) on the lines of , with respect to the scalar product [START_REF] Gross | Two methods for Wavelet-Based volume rendering[END_REF], in order to get orthonormal basis functions. Gram-schmidt orthogonalization process is known to be numerically instable for big matrices, but our matrices are always of small dimension (the dimension of our matrices is bounded by the degree of the removed vertex). Moreover orthonormality is crucial to ensure the stability of A u t h o r p r e p r i n t the algorithms: since the coarsest coefficients can result from millions of analysis matrices multiplications of the finest coefficients for big data sets, using non-orthonormal analysis matrices leads to instable algorithms.

To sum up, the computation of the analysis matrices and requires the following steps:

Basic hierarchical decomposition

As previously stated, the local decomposition described in section 2 can be applied together with various triangle decimation algorithms.

In this section we apply this local decomposition together with a hierarchical structure used in several previous works: the hierarchical Delaunay triangulation.

We begin by some background on hierarchical Delaunay triangulation and give the reasons for this choice. Then we show two series of examples that aim to test the stability of our method, as well as the validity of our choice for the wavelet components.

Hierarchical Delaunay triangulation.

Hierarchical Delaunay triangulations were heavily used in previous works, especially in the field of computational geometry. They were introduced in 1983 by Kirkpatrick [START_REF] Kirkpatrick | Optimal search in planar subdivisions[END_REF] in the purpose of point location algorithms. They were later successfully applied in the field of terrain modeling ( [START_REF] De | A pyramidal data structure for trianglebased surface description[END_REF]).

Starting from a surface triangulation, the construction of the hierarchy can be sketched as follows:

While (# selected vertices ) select a set of at least independent (non-adjacent) vertices with degree less than (where is the current number of vertices, and and are constants: and were used in Kirkpatrick's paper), remove each of the selected vertex and their adjacent edges creation of blocks, retriangulate the blocks according to a Delaunay criterion, update the links between the levels. Kirkpatrick proved that the hierarchy can be built in linear time, and has a logarithmic depth. This algorithm was introduced for planar triangulations, but can be applied on more general surface geodesic triangulations, with the assumption that star-shaped polygons on the surface can be triangulated. This is in particular the case for all convex surfaces. We will show planar and spherical triangulations. Our current implementation is based on a greedy algorithm in a first step, followed by a local swap of the diagonals in each quadrilateral, in order to ensure the convexity of the triangulation for spherical triangulations, or an empty circle property for planar triangulations. This algorithm can be implemented with several data structures. Depending on the data structure, the links can join either faces, vertices or blocks between different levels. Fig. 6 shows an example of such a hierarchy. The original and simplified triangulations are shown on the left. Empty circles point to vertices that are removed to reach the next level. The middle part shows the corresponding hierarchies over the faces. Black circles denote triangles, and links join triangles that have non-empty intersections. Roundedrectangles group triangles belonging to the same block. The right part of fig. 6 shows these blocks, and the links between intersecting blocks. Each block has two triangulations: triangles before the vertex removal (drawn solid) are called child triangles of the block, and triangles after the vertex removal (in dashed) are called parent triangles. Since only one vertex is removed in each block, the number of wavelet coefficients generated in each block equals two (difference between the number of parent and child triangles in each block). Moreover, these two wavelet coefficients can be assigned to the removed vertex. This last feature will be useful in section 4, in order to improve the approximations.

After the whole analysis of the data set, each removed vertex has been assigned two wavelet coefficients, and the resulting data set is the coarsest approximation of the original one.

Level & threshold reconstructions

Once the hierarchy is constructed, and the wavelet coefficients for each block are computed, then two different kinds of partial reconstruction may be performed.

Level reconstruction

Level reconstructions result from the application of the reconstruction formula (4) on each block of the hierarchy, starting at the root level, and going down to the desired final level of reconstruction. Level reconstructions do not need the knowledge of the links between blocks. The only required information is the list of blocks, ordered from the root level to the finest level.

Threshold reconstruction

The insertion of one wavelet coefficient in the approximation implies obviously the insertion of the vertices of the block from which that wavelet coefficient has been generated. But it requires also the insertion of the blocks in the upper levels of the hierarchy, which have a link (e.g. a non-empty intersection) with that block. The necessity of including triangles that lie outside the support of the wavelet is not a consequence of using hierarchical Delaunay triangulations. In fact, the same result holds for 4-to-1 split triangulations. Fig. 7 illustrates this necessity. The top part shows the hierarchies. The solid rounded-rectangle denotes the block from which a wavelet coefficient has to be inserted. The dashed roundedrectangles show the blocks that have to be inserted because of the hierarchy, in order to avoid hanging nodes. The bottom part of fig. 7 shows the corresponding triangulation: Solid triangles are the ones which triangulate the support of the wavelet. Dashed triangles are required in order to avoid hanging nodes.

Therefore we see that the partial reconstruction of the data set, using wavelet coefficients greater than a given threshold can be performed with the following steps:

select the minimal sub-hierarchy that contains the wavelet coefficients greater than the given threshold insert the blocks in the sub-hierarchy, starting from the root level, and going down to the finest level of the selected wavelet functions.

Figure 7: Including a wavelet coefficient: triangles lying outside the wavelet support are required.

Numerical results

In this section we show some numerical examples of analysis/reconstruction on a data set that was specifically constructed to test the stability of the algorithms. We have chosen a strongly irregular planar triangulation, and have mapped a non-smooth image on that triangulation: the triangulation was generated by Delaunayinserting random points in a square domain, and the data value on each triangle equals the mean value of an image generated from the ETOPO5 data set from the National Oceanographic and Atmospheric Administration 1 , showing the elevation (depth) of the earth.

The data set consists of 50002 faces. The resulting hierarchy has 32 levels. Fig. 8 shows the logarithmic plot of the number of triangles versus level index. As expected, the curve is linear, illustrating the results of Kirkpatrick ([9]). Fig. 9 shows the relativeerror versus number of triangles for the two partial reconstructions described in section 2.2. The approximations are clearly better by threshold reconstruction than by level reconstructions, thus illustrating the validity of our choice for the wavelet coefficients. The effect of choosing the biggest wavelet coefficients is also illustrated on fig.

Improving the approximations

In this section we investigate two ways of improving the approximations resulting from the decomposition/reconstruction algorithms described in section 3. The first idea is to set some constraints during the construction of the pyramidal structure, in order to preserve important features of the data. And the second idea is to use the wavelet coefficients to construct a second pyramid that gives better approximations than the initial one. These ideas are illustrated on "real size" spherical data sets at the end of this section. 1 Available via anonymousftp at ftp://ftp.ngdc.noaa.gov/Solid Earth/Top A u t h o r p r e p r i n t

Preserving important features

Using the -norm as a distance error may sometimes result in the lost of important features in the data set, even for an approximation very close to the original data. Think of the coastline in an elevation/depth data set for the earth: if the coast is sharp it will be preserved in a close -approximation, but if it is smooth (between England and France for example), it will be moved, or it can even disappear.

A simple yet efficient way to preserve such important features is to mark vertices, in order to forbid their removal during the construction of the hierarchal structure. We have done this for the earth elevation/depth data set (see fig. 11): all vertices along the coastlines at the finest level have been preserved during the analysis. This implies that all wavelet functions have a support that doesn't intersect the coastlines. Therefore, reconstruction with any threshold exactly reproduces the data set along the coastlines. This approach is quite different from the one described in [START_REF] Schroeder | Spherical wavelets: texture processing[END_REF]: since their wavelet supports, defined on subdivision meshes, can intersect the coastlines, they have to check all such wavelet functions, and to require the insertion of all corresponding detail coefficients.

Our method has been used also in order to handle the data set in color plate 2: as explained in detail in section 4.3, this data set is defined on a restricted spherical domain. We were able to handle this data set by forbidding the removal of all vertices outside, and along the domain of the data set.

Optimizing the hierarchical structure

We already have pointed out in section 3.2 that the insertion of one wavelet coefficient implies the insertion of triangles that lie outside the support of the wavelet, in order to avoid hanging nodes. These unwanted triangles are in the upper (coarser) levels of the hierarchy, and belong to blocks that have a non-empty intersection with the wavelet support. Therefore the ideal hierarchical structure would have biggest wavelet coefficients in the upper (coarser) levels of the hierarchy, and smallest coefficients in the lower (finer) levels. Such a hierarchical structure would give optimal results, in the sense that the insertion of wavelet coefficients in decreasing order of magnitude would not require the insertion of unwanted triangles. However, constructing such an optimal hierarchy would be computationally too expensive. In order to still improve the approximations for a reasonable computational cost, we have done the following: compute the two wavelet coefficients for each vertex as explained in section 3, sort the list of removed vertices according to the sum of their squared wavelet coefficients (lowest first, biggest last), construct a new sequence of approximations from the original data, and removing each vertex in the order given by the previ-A u t h o r p r e p r i n t ous step. Thereby two new wavelet coefficients are computed for each removed vertex.

This preprocess requires two full analysis of the data set instead of one, as well as sorting the list of vertices. On the other hand, the resulting improved sequence of approximations is easier to reconstruct: the vertices are inserted one after another together with their wavelet coefficients, regardless of the norm of these wavelet coefficients. No hierarchy has to be traversed in order to avoid hanging nodes.

Fig. 12 compares the threshold reconstruction in the numerical example of section 3.3 with the new sequence of approximations. Fig. 13 shows the resulting new approximation with the same number of triangles as in fig. 10. Note the improved concentration of the vertices in discontinuous areas of the data set.

Results

Implementation issues

In the analysis of space and time requirements of our algorithms, it is important to distinguish the geometrical part from the numerical part: the geometrical part includes computing and storing the Delaunay hierarchy of meshes, whereas the numerical part involves computing and storing the coarse and detail coefficients.

As stated in section 3.1, the geometrical part has been extensively worked in others papers. [START_REF] De Berg And Katrin | On levels of detail in terrains[END_REF] for ex. describes a data structure that is based on blocks, as defined in section 3.1. The size of all proposed data structures is proportional to the number of vertices, with a bigger constant if the adjacencies are stored. The computation time for the construction of the hierarchy (without analysis of the data) is also proportional to the number of vertices. The highest cost operation is the local retriangulation of the blocks. Our current implementation, which is not optimized, works at about 3000 vertices per second on a SGI Octane 175Mhz R10000.

The numerical part of our algorithm requires the storage of the coarse and detail coefficients. As explained in section 2.2.1, this can be done with a constant memory size, since these coefficients replace the fine coefficients at each local analysis. The analysis and synthesis matrices are not stored in memory. The computation time for the numerical part is mainly due to the evaluation of the intersections between the coarse and fine triangles. Again this has not been optimized in our current code, and it works at about 4000 vertices per second on a SGI Octane 175Mhz R10000.

Examples

Color plate 1 illustrates sections 4.1 and 4.2 on a spherical data set consisting of 1.3M faces. The original data set is defined on a regular triangulation (4-to-1 split with 8 levels of subdivision starting on a icosahedron). The data value on each face equals the mean value of the ETOPO5 data set (which consists of 2160 x 4320 samples on a uniform grid). A first full analysis has been performed on the data set, with preservation of the coastlines, as explained in section 4.1. (a), (b) and (c) show three snapshots taken during this first analysis. (c) shows the coarsest level: no more vertices can be removed, the only remaining vertices are those along the coastlines, and those with degree greater than 11. Although the faces are rendered as planar, all computations are done with true spherical triangles (in particular areas and intersection computations). The vertices were sorted according to the wavelet coefficients, and a second full analysis was performed, as explained in section 4.2. (d), (e) and (f) show three snapshots taken during the reconstruction.

Color plate 2 shows the gravity anomaly on the earth. The gravity anomaly is the difference between the actual gravity, and the ideal gravity for a perfect ellipsoidal surface. The data set comes from satellite measurements. For technical reasons, it is restricted to the water surface on the earth, and between and degrees of latitude. More information on this data set, as well as the data itself are available at http://topex.ucsd.edu/marine grav/mar grav.html. The original data consists of short numbers. Even numbers signify no measurement available while odd numbers contains the gravity anomaly . Our method based on irregular grids can handle such complicated domains. Starting with the same regular grid as for color plate 1, we have used the method described in section 3.1 in order to remove almost all vertices on the continents and below or above degrees of latitude. Of course we did not store the hierarchy for this initial sequence of vertex removals. The resulting irregular triangular mesh consists of 920250 triangles, with small triangles covering exactly the data domain, and some few big triangles covering the areas where no data is available. Our algorithm has been applied on this initial triangulation, using the method described in section 4.1 in order to preserve the domain boundaries (i.e. coastlines and parallels).

Conclusion

We have introduced algorithms for decomposing and reconstructing piecewise constant functions defined on irregular triangular meshes. Our algorithms combine hierarchical Delaunay triangulation with a wavelet-like decomposition in order to build a multiresolution model of the data set. Our local decomposition generalizes Haarwavelets for non-nested triangulations, and is therefore especially adapted for the analysis of big and irregular data sets. Future work will address the generalization to higher order wavelets, in order to handle smoother data sets.

Figure 1 :

 1 Figure 1: Local triangle decimation.

Figs. 2

 2 Figs. 2 and 3 illustrate our local decomposition on two examples. In both figures, the top part shows from left to right, the finer function , the intermediate function , and the final coarser function . The bottom part shows the detail coefficients times the wavelet functions: . Fig. 2 shows the local decomposition on a 4-to-1 split example. This leads to a traditional Haar wavelet decomposition. Note in this case that the intermediate function (top-middle),

Figure 2 :

 2 Figure 2: Local decomposition by 4-to-1 split: finer, intermediate and coarser approximations on top, detail coefficients times wavelet functions on bottom. In this case, the coarse and fine triangular domains are nested, and therefore, the intermediate approximation equals the coarser approximation. The relative high magnitudes of the detail coefficients (bottom part) show the big error between the fine and coarse approximations.

Figure 3 :

 3 Figure 3: Local decomposition by Delaunay-removal: finer, intermediate and coarser approximations on top, detail coefficients times wavelet functions on bottom. Since the triangular domains are nonnested, the intermediate approximation differs from the coarse approximation. The relative low magnitude of the detail coefficients (bottom part) show the small error between the fine and coarse approximations.

Figure 4 :

 4 Figure 4: Computation of : intersections between finer and coarser triangles.

  computation of the intersections between coarse and fine triangles, computation of the relative areas area area orthogonalization process on a matrix.

Fig. 5

 5 Fig. 5 illustrates one step on the construction on a planar example: empty circles on the left denote vertices that are selected for removal. The right part shows the result after retriangulation.

Figure 5 :

 5 Figure 5: One step of the construction of the hierarchy: (a) independent vertices selection, (b) retriangulation.

Figure 6 :

 6 Figure 6: Hierarchical Delaunay triangulation

10 :

 10 two snapshots were taken, with almost the same number of triangles. The top part shows the reconstruction at level 21, consisting of 2326 triangles (top left=data set, top right=triangulation). The bottom part shows a threshold reconstruction consisting of 2388 triangles (bottom left=data set, bottom right=triangulation). Note how the selection of the biggest wavelet coefficients affects the triangulation, adding more triangles in discontinuous areas of the data set.

Figure 8 :

 8 Figure 8: level vs #triangles.

Figure 9 :

 9 Figure 9: Relative -error vs. number of triangles. 'o': level reconstruction, '+': threshold reconstruction.

Figure 10 :

 10 Figure 10: Top=level 21 reconstruction, 2326 triangles (left=data set, right=triangulation), Bottom=threshold reconstruction, 2388 triangles (left=data set, right=triangulation).

Figure 11 :

 11 Figure 11: A coarse reconstruction (120000 faces out of 1.3M), with (left) and without (right) preservation of the coastlines. Although the error is smaller without preservation of the coastlines, the result is clearly not acceptable for this data set.

Figure 12 :Figure 13 :

 1213 Figure 12: Relative error vs. number of triangles: '+'=threshold reconstruction with non-optimized pyramid, 'x'=reconstruction with improved sequence of approximations (see 4.2)

  (a), (b), (c) and (d) show three snapshots taken during the reconstruction with the optimized pyramid. Deep blue indicates low gravity values, while red indicates high gravity values. (e), (f), (g) and (h) are close-up views above the pacific A u t h o r p r e p r i n t ocean, where the gravity anomaly is especially irregular ((e) and (h) show the edges of the triangulation in light blue).

Acknowledgments

The author wishes to thank the referees for their helpful comments.


Georges-Pierre Bonneau

Alexandre Gerussi

Level of detail visualization of scalar data sets on irregular surface meshes

Keywords: wavelets, non-regular triangulations, compression, visualization

, where a similar framework is discussed for piecewise constant functions.

INTRODUCTION

Our aim is to build a sequence of level of detail representations of continuous piecewise linear scalar data sets defined on irregular planar or spherical triangulations. Our approach is based on vertex removal operations, and combines a local decomposition/reconstruction algorithm inspired from wavelet techniques, together with a global greedy algorithm, in order to construct the levels of detail. Following our previous work dedicated to piecewise constant data sets ( [START_REF] Bonneau | Hierarchical decomposition of datasets on irregular surface meshes[END_REF]), we use a local wavelet-like decomposition, that maps a functional space onto a coarser space. The fact that the coarser space has no to be a subspace of the finer one (which must be the case in wavelet theory) enables to deal with irregular triangulations. At the global level, the choice of the order of removal of the vertices is made by a greedy algorithm that is guided by a simple error criterion.

In section 2, we briefly explain the global greedy algorithm. Section 3 describes the functional spaces and the basis functions result-LMC-CNRS, BP.53, F-38041 Grenoble Cedex 9, France.

ing from the vertex removal process. Section 4 is dedicated to the local decomposition/reconstruction at each vertex removal. Eventually section 5 presents numerical and visual results on a planar data set, and on a spherical data set with 1.3M faces.

Previous work.

Previous work on level of detail representations are mostly dedicated to surface simplifications, i.e. they deal with geometric data sets, while we restrict ourselves with scalar data sets ( [START_REF] Hoppe | Progressive meshes[END_REF][START_REF] Hoppe | View-dependent refinement of progressive meshes[END_REF][START_REF] Garland | Surface simplification using quadric error metrics[END_REF]). In other words, we are not trying to approximate a surface, we approximate a scalar function defined on a fixed surface.

Schröder and Sweldens ( [START_REF] Schröder | Spherical wavelets: Efficiently representing functions on the sphere[END_REF]) have introduced spherical wavelets to build sequences of level of detail approximations to scalar data sets on the sphere. But since their approach relies strictly on the framework of wavelet theory, it is restricted to data sets defined on regular triangulations constructed by recursive 4-split of a base mesh.

Papers on level of detail for TIN (Triangular Irregular Network) deal with scalar data sets (height values) on irregular planar triangulations (see for example [START_REF] De | A pyramidal data structure for triangle-based surface description[END_REF]). [START_REF] Chandrajit | Error-bounded reduction of triangle meshes with multivariate data[END_REF] deals with multivariate scalar data sets on arbitrary surface meshes. Like our paper, they are mostly based on vertex removal. But the crucial difference is that, after each vertex removal, we modify the values of the adjacent vertices, while they leave these values unchanged. In section 5, we will give some visual and numerical experiments that compare sub-sampling results with our algorithm, for the same sequence of vertex removal.

Global level algorithm

Our algorithm is based on vertex removal operations. At the global level, a greedy algorithm based on a simple cost function is used to decide the order in which the vertices are removed. If M is a vertex with scalar value , and Mi are the vertices adjacent to M, with value i, we use the following cost for the removal of vertex M:

cost of removal = P i 1 + distmax ; d(M Mi) i P j 1 + distmax ; d(M Mj) ; (1)
where d(M Mi) is the geodesic distance between M and Mi, and distmax = m a x i d(M Mi).

We could have used the exact L2 error induced by the removal of each vertex, since our local decomposition is based on L2 approximation. But we found that the simple cost function (1) reduces the computation time noticeably and still works well in practice.

In addition to the cost function [START_REF] Chandrajit | Error-bounded reduction of triangle meshes with multivariate data[END_REF], it is possible to forbid the removal of some vertices in order to preserve important features in the data set. We have used this technique to preserve coast-lines in the earth data set (see section 5).

A u t h o r p r e p r i n t 3 Functional spaces and basis functions

As stated in the introduction, we want to deal with piecewise linear functions. In fact in this paper we'll show results for such functions in two cases: a planar triangular mesh and a spherical one. The spherical case can be treated almost as the planar one, except for some details (for example the distances on the sphere are the geodesic distances, the re-triangulation criteria is based on a convex-hull property). We will now implicitly consider, in our illustrations, planar meshes.

We start with a 2D-triangulation Tn of n vertices M1 : : : M n. The global greedy algorithm described in section 2 is used to compute a sequence Tn T n;1 : : : T m of simplified triangulations. The domain covered by Ti will be denoted by D(Ti). By removing a vertex M from the mesh, the triangulation is only modified locally, over the so-called polygon of influence (PI) of M, which is the polygon whose vertices are adjacent to M (the 1neighbours of M). This is illustrated in figure 1.

M PI(M)

Figure 1: removing of a vertex from a triangulation

There are several ways to choose the new edges which re-triangulate the polygon. Our current implementation is based on a local greedy algorithm in a first step, followed by a local swap of the diagonals in each quadrilateral, in order to ensure the convexity of the triangulation for spherical triangulations, or an empty circle property for planar triangulations.

The sequence Ti, i = m : : : nis used to define the functional spaces forming our multi-resolution analysis: Let Vi, i = m : : : n , be the space of all functions that are defined and continuous over D(Ti) and piecewise linear over each triangle of the mesh Ti (shortly: a CPLF, for continuous piecewise linear function).

The definition of a function of Vi is equivalent to the knowledge of its values at every vertex of Ti, and consequently the dimension of Vi is i. The basis of Vi is given by the i hat-functions: the hatfunction of the vertex Mj in Vi takes value 1 at Mj and then declines linearily to zero on each triangle of the PI of M in Ti (see figure 2 for an illustration).

Local decomposition/reconstruction

The problem of computing the decomposition/reconstructionis similar at each level. Consequently, we shall now focus on the transition between Vi and Vi;1 for a fixed i, or, to hide the index i, between V f and V c (f=fine, c=coarse), where V c results from the removal of vertex M in V f . Our decomposition/reconstruction process is based on L2approximation between V f and V c . More precisely, given a function F f 2 V f , we want to compute two things: This decomposition/reconstruction process is inspired from wavelet-based techniques, and extents them, in the sense that it can be applied even if the coarse space V c is not included in the fine space V f . This extension has been introduced in a previous work in a different context ([3], section 2). The reader might refer to [START_REF] Bonneau | Hierarchical decomposition of datasets on irregular surface meshes[END_REF] or [START_REF] Bonneau | BLaC-Wavelets: A multiresolution analysis with nonnested spaces[END_REF] for a context-independent presentation of this extension.

its L2-approximation F c in V c (which is the orthogonal pro- jection of F f in V c ),

Basis functions in V f and V c

Before going further with the actual computations, we should have a careful look at the basis functions in V f and V c . The hat function corresponding to the vertex M no longer exists in V c . This is the reason why dim(V f ) = dim(V c ) + 1 : Also, the hat functions based on the neighbours of M in T f disappear, since the vertex removal creates a new neighbourhood for all of them, but they are "replaced" by hat functions defined on the same vertices and their new neighbourhoods in T c . All the other basis function remain unchanged. By comparing the supports of the basis functions in figure 3, you can see the consequences of a vertex removal on the basis functions whose support overlaps the PI of the removed vertex.

Before vertex removal

After vertex removal We note M f i , i = 0 : : : kand M c i , i = 1 :: k the basis of V f and V c , respectively. M f i is the hat function centered on vertex Mi in T f . We have just seen that, in fact, if we number 1 to l the 1neighbours of M (in T f ) then M f i = M c i 8i = l + 1 : : : k :

Unchanged basis function

Modified basis function

A u t h o r p r e p r i n t

To perform L2-approximation between V f and V c , we will have to compute the following matrices of scalar products between finer and coarser basis functions:

G f = ( < M f i M f j >)i=0 ::: k j=0 ::: k G c = ( < M c i M c j >)i=1 ::: k j=1 ::: k U cf = ( < M c i M f j >)i=1 ::: k j=0 ::: k A = ( G c ) ;1 U cf G c and G f are the Gram-schmidt matrices in V c and V f respectively. And for those who are familiar with the terminology of wavelet theory, A can be seen as the matrix of scalar products between the fine scaling functions (in V f ) and the dual coarse scaling functions (in V c ).

Decomposition

The decomposition step takes as input a function

F f = P k i=0 a f i M f
i in the finer space V f , and outputs an approximation F c = P k i=0 a c i M c i in V c and a detail coefficient b. Note that there is only one detail coefficient since dim(V f ) = dim(V c ) + 1 . This detail coefficient is conceptually linked to the vertex M, it can be stored in place of a f 0 .

Coarse coefficients

The coarse coefficients (a c : ) are computed from the fine coefficients (a f : ) with the following formula:

(a c : ) = A(a f : ) :

In fact the matrix A ensures that F c is the best L2-approximation of F f in V c .

Detail coefficient

The detail coefficient is computed from the (a f : ) coefficients through the equation: b = B(a f : ) where B is a line-vector representing a dual wavelet function, in the basis of V f (see [START_REF] Bonneau | Hierarchical decomposition of datasets on irregular surface meshes[END_REF] or [START_REF] Bonneau | BLaC-Wavelets: A multiresolution analysis with nonnested spaces[END_REF]). This vector is computed as follows: First we know that it has to fulfill the orthogonal property A ; G f ;1 B T = 0 (which means in analogy with wavelet theory, that the dual wavelet is orthogonal to the dual coarse scaling functions). This gives a one-dimensional space of solutions for the line-vector B. B is then determined uniquely by norming it to 1 in V f : B ; G f ;1 B T = 1 . The orthogonalization and normalization on B ensure that the detail coefficient b is a measure of the approximation error between the input function F f and the output function F c .

Localization

One problem of importance is that the matrices G c , G f and U cf have to be computed at each vertex removal, and this is a too big amount of work since these matrices, although sparse, can potentially be very large. Nevertheless, since a vertex removal only changes basis functions in a small neighbourhood, it is quite intuitive that the best approximation of F f is not going to differ very much from F f outside of this neighbourhood.

Precisely, we parameter our decomposition/reconstruction process with a integer K 1 that reduces the global problem to a local one by doing the following: Instead of working with all basis functions, we only consider the basis functions which are based on vertices which are i;neighbours of M in T f , i = 1 : : : K(see figure 4). This means that F c is no longer the best L 2 -approximation of F f but it is supposed to be very close to. We won't try to quantify this difference here, but this assumption led to good results (K = 2 or 3 appeared to be large enough).

1-neighbourhood

2-neighbourhood

Figure 4: the 1-and 2-neighbourhood of a vertex Another consequence of this strategy is that the global algorithm is now in o(n), because each local analysis/synthesis is made in constant time, independently of the number of vertices in Ti.

Computing the scalar products

Let's say a few words about the computations of the matrices G c , G f and U cf . Their coefficients are scalar products between the functions M f : and M c : . Because M f i = M c i whenever i > l , these three matrices have a lot of coefficients in common. Moreover, a function only has non zero scalar product with its 1-neighbour functions, and the support of two such functions only overlap over two triangles (see figure 5).

Figure 5: intersection of the supports of two basis functions.

In fact, the only difficult cases occur in the computation of scalar products of mixed type < M f i M c j > when i l and j l, because it is then necessary to compute the intersections between the fine and coarse triangles forming the inside of the PI of M in T f and T c respectively. Nevertheless, without sinking into too much detail, this problem is simpler than the general problem of finding the intersections between any two triangle sets, because the PI is a star-polygon, centered on the removed vertex.

Finally, the computation of the matrices only requires a reasonable amount of work, and this is crucial because our framework is intended to run on possibly very large datasets (see the next section and the application on the 1.3M faces spherical data set).

Reconstruction

The reconstruction process is relatively straightforward: the synthesis matrix is nothing but the inverse of the analysis matrix, thus, the A u t h o r p r e p r i n t fine coefficients (a f : ) can be recovered from the coarse coefficients (a c : ) and the detail coefficient b with the following formula:

(a f : ) = A B ;1 a c : b

Illustration

We shall conclude section 4 with an illustration of our decomposition algorithm. Figure 6 shows at the top left an input function (the removed vertex is circled), at the top right the output function, and the bottom shows two views of the corresponding wavelet function. This decomposition was computed using the localization parameter

K = 2 .

Results

We start this section with some numerical examples of decomposition/reconstruction on a data set that was specifically constructed to test the stability of the algorithms. We have chosen a strongly irregular planar triangulation, and have mapped a non-smooth image on that triangulation: the triangulation was generated by Delaunayinserting random points in a square domain, and the data value at each vertex was sampled from the ETOPO5 data set from the National Oceanographic and Atmospheric Administration 1 , showing the elevation/bathymetry of the earth.

The data set has 25000 vertices. The coast-lines were preserved by the global greedy algorithm, as explained in section 2: the vertices to keep were detected by looking if a change of sign occured at the 1-neighbour vertices. The decomposition was computed using the localization parameter K = 2 (see section 4.3). We compare the results of our decomposition algorithm with the results obtained by simply sub-sampling the data set. Figure 7 shows the relative L 2error versus number of vertices for our decomposition, and for subsampling. Figure 8 shows two partial reconstructions with 2500 vertices out of 25000, one based on our algorithm (bottom right), and the other based on sub-sampling (bottom left), for the same sequence of removed vertices.

Color plate 1 illustrates our algorithm on a spherical data set consisting of 1.3M faces. The original data set is defined on a regular triangulation (4-to-1 split with 8 levels of subdivision starting on a 1 Available via anonymousftp at ftp://ftp.ngdc.noaa.gov/Solid Earth/Topography/tbase 5min/ icosahedron). The data value at each vertex was sampled from the ETOPO5 data set (which consists of 2160 x 4320 samples on a uniform grid).

Our algorithm was applied using the localization parameter K = 2

Figure 2 :

 2 Figure 2: Example of a CPLF and its decomposition in the hat function basis

Figure 3 :

 3 Figure 3: Supports of hat basis functions, before and after vertex removal.

Figure 6 :

 6 Figure 6: Approximating a CPLF: top left = input function, top right = approximation, bottom = wavelet function (two views). The removed vertex is circled.

(see section 4 . 3 )

 43 , and with preservation of the coast-lines. The upper part of color plate 1 compares the result our algorithm, with the result obtained from sub-sampling the data set: (a) shows the approximation resulting from our algorithm, and (c) shows the sub-sampling result. (b) shows the corresponding spherical mesh. The bottom part of color plate 1 shows different partial reconstructions computed from our algorithm: (g), (h), (i) show the spherical meshes corresponding to (d) (100000 vertices), (e) (200000 vertices) and (f) (300000 vertices).

Figure 7 :Figure 8 :

 78 Figure 7: Relative L2-error vs. # of vertices. 'o':subsampling,'+':L2-approximation


Georges-Pierre Bonneau
email: [bonneau@imag.fr

Alexandre Gerussi
email: alexandre.gerussi]@imag.fr

Hierarchical decomposition of datasets on irregular surface meshes

In this paper we introduce multiresolution analysis (MRA) algorithms intended to be used in scientific visualization, and based on a non-nested set of approximating spaces. The need for non nested spaces arises from the fact that the required scaling functions do not fulfill any refinement equation. Therefore we introduce in the first part the concept of approximated refinement equation, that allows to generalize the filter bank and exact reconstruction algorithms. The second part shows how this concept enables to define a MRA scheme for piecewise constant data defined on an arbitrary planar or spherical triangular mesh. The ability to deal with arbitrary triangular meshes, without subdivision connectivity, can be achieved only through the use of non nested approximating spaces, as introduced in the first part.

Introduction

Hierarchical decomposition of datasets based on wavelet theory is a growing subject of research in computer graphics and scientific visualization. Wavelet theory is based on the knowledge of a sequence of functional spaces in which the data is successively approximated. Unfortunately, this sequence has to be nested. This explains why wavelet theory can't be applied to the multiresolution analysis of datasets defined on irregular meshes, since such meshes cannot be reach by a subdivision process starting on a coarser mesh. Recently ( [START_REF] Bonneau | BLaC-Wavelets: A multiresolution analysis with non-nested spaces[END_REF]) we have introduced a theoretical framework that enables to deal with a non-nested sequence of functional spaces. The link between successive non-nested spaces is given through the definition of intermediate spaces. We will show how this framework can be applied to the multiresolution analysis of datasets defined on irregular surface meshes. We will combine our results with a well-known hierarchical structure-the hierarchical Delaunay triangulation-in order to obtain successive approximations of the original dataset on simplified meshes. We will show several examples that demonstrate the stability of our decomposi-tion/reconstruction algorithms, as well as the validity of our choice for the wavelet functions. We will explain how we can use the wavelet coefficients in order to optimize the hierarchical decomposition of the dataset.

Multiresolution Analysis with non-nested approximating spaces

We assume that the reader is familiar with the basic concepts of multiresolution analysis. An introduction on this topic can be found in [START_REF] Sollnitz | Wavelets for Computer Graphics: Theory and Applications[END_REF]. In this section, a general framework for constructing multiresolution analysis schemes based on a non nested sequence of approximating spaces is given. All functional spaces throughout the paper are supposed to have finite dimensions.

Suppose you are given a sequence of functional spaces V n that satisfies the following condition:

V n is isomorph to a subspace of V n+1 : [START_REF] Bonneau | BLaC-Wavelets: A multiresolution analysis with non-nested spaces[END_REF] This condition means that the sequence of approximating spaces is "growing": it may be rewritten as dim V n dim V n+1 .

Approximated refinement equation

Let (' n j ) denote the scaling functions (i.e. (' n j ) is a basis of V n ). If the usual nested condition V n V n+1 would be satisfied, the link between two successive levels of resolution would be provided by a so-called refinement equation, which expresses each scaling function ' n j as a linear combination of the finer scaling functions (' n+1 i ). In our setting, condition (1) implies that there exists functions ( ' n j )

which form a basis of a space Ṽn , isomorph to V n , and included in V n+1 . And the link between level n and n + 1

is now provided through a two step process: First apply the isomorphism, i.e. replace the scaling functions (' n j ) by the intermediate functions ( 'n j ), and then apply the refinement equation which expresses 'n j as a linear combination of the finer scaling functions (' n+1 i ):

' n j $ 'n j = X i p n ij ' n+1 i :

(

We refer to [START_REF] Kirkpatrick | Optimal search in planar subdivisions[END_REF] as the approximated refinement equation A crucial point of our method is the choice of the intermediate scaling functions ( 'n j ). This choice will be precised later.

Wavelet functions and filter bank

Once the intermediate functions, and thus the intermediate spaces Ṽn are chosen, we choose the detail space W n as a complementary of the intermediate space Ṽn in the finer space V n+1 : V n+1 = Ṽn W n :

The wavelet functions ( n j ) are basis functions of the space W n . Thus there exist coefficients (q n ij ) such that

n j = X i q n ij ' n+1 i and coefficients (a n ij ), (b n ij ) such that ' n+1 j = X i a n ij 'n i + X i b n ij n i :
Let A B P Q denote the matrices with coefficients

(a n ij ) (b n ij ) (p n ij ) (q n ij ).
Given an approximation

f n+1 = P i x n+1 i ' n+1 i in V n+1
, the next coarser approximation f n = P i x n i ' n i and the details g n = P i y n i n i are computed by

(x n i ) = A(x n+1 i ) (y n i ) = B(x n+1 i ) : (3)
The inverse transformation is given by

(x n+1 i ) = P(x n i ) + Q(y n i ) : (4)
The analysis and reconstruction formulas (3) and ( 4) are the same as if the nested condition V n V n+1 would be satisfied, but they are related now to an approximation in two steps:

f n+1 app. ! e f n = P i x n i e ' n i iso. $ f n = P i x n i ' n i & g n

Choice of the analysis and synthesis matrices

The complete choice for the analysis and synthesis matrices depends on the application, and thus will be made in the next section. Nevertheless, some conditions that these matrices have to fullfill in order to ensure best-approximation and orthogonality properties can be precised now.

Let G n denote the Gram-schmidt matrix of the basis ' n i (i.e. the matrix whose elements are < ' n i ' n j >). If we want f n , the approximation at resolution n, to be the best L 2 approximation of f n+1 , then it is easy to show that the analysis matrix A must satisfy the following condition:

A = G ;1 n (< ' n i ' n+1 j >) : (5)
Once A is computed from ( 5), the matrix B can be chosen orthogonal to A, and orthonormal, with respect to the scalar product in the dual scaling function basis of V n+1 :

AG ;1 n+1 B T = 0 (6)
BG ;1 n+1 B T = I :

P and Q are then computed from A and B in order to ensure the exact reconstruction property:

( P Q ) = A B ;1 : (8)
Together, conditions ( 5), ( 6), (7), and (8) ensure that the approximation at resolution n f n is the best L 2 approximation of both f n+1 and the intermediate approximation fn .

Application to piecewise constant data sets on irregular triangulations

In this section, the theoretical results of section 2 are applied to the multiresolution analysis of piecewise constant data sets on irregular planar or spherical triangular meshes.

Previous works on wavelet methods for data sets defined on triangular meshes are restricted to regular meshes obtained by recursively spliting in 4 triangles each triangle of a coarse base mesh ( [START_REF] Schröder | Spherical wavelets: Efficiently representing functions on the sphere[END_REF], [START_REF] Nielson | Haar-wavelets over triangular domains with applications to multiresolution models for flow over a sphere[END_REF], [START_REF] Lounsbery | Multiresolution analysis for surfaces of arbitrary topological type[END_REF]) (see fig. 1). This recursive process is associated to a simple hierarchy, that has the structure of a forest of quadtrees. For such meshes, it is possible to construct a nested sequence of approximation spaces.

Figure 1. 4-to-1 split

When dealing with irregular triangulations, it is not possible to associate a simple tree structure to the data set. Other possible hierarchies may be constructed by applying edge contraction, triangle contraction or vertex removal operation on the triangular mesh. These hierarchies cannot be represented by tree structures, and no sequence of nested approximation spaces can be associated. We have chosen to use vertex removal operations (see fig. 2), and more specifically we have applied previous works from Kirkpatrick ([2]) to build a so-called hierarchical Delaunay structure. This is explained in subsection 3.1.

Figure 2. Point removal

Once the hierarchy has been computed, we have applied the results of section 2 locally to each star-shaped polygon created by the vertex removal operations, in order to compute new face values on the new triangles, as well as detail values that measure the distance error between the original data set and the new simplified one (as illustrated in fig. 3). This will be explained in subsection 3.2. y 0 y 1

Figure 3. Computation of new face values and detail coefficients

Eventually we show in subsection 3.3 how the local analysis and synthesis described in 3.2 combine with the hierarchy introduced in 3.1 in order to yield global analysis and synthesis algorithms for piecewise constant data sets defined on irregular planar or spherical triangular meshes.

Hierarchical Delaunay T riangulation

Starting from a surface triangulation, the construction of a hierarchical Delaunay structure, as introduced in 1983 by Kirkpatrick ([2]), can be sketched as follows: While (# selected vertices > 0) select a maximal set of independent (non-adjacent) vertices remove each of the selected vertex and their adjacent edges ) creation of blocs, retriangulate the blocs according to a Delaunay criteria, update the links between the levels.

Although it was introduced for planar triangulations, this algorithm can be applied to spherical triangulations as well. It can be implemented with several data structures. Depending on the data structure, the links can join either faces, vertices or blocs between different levels. Figure 4 shows an example of such a hierarchy. The original and simplified triangulations are shown on the left. Empty circles point to vertices that are removed to reach the next level. The middle part shows the corresponding hierarchies over the faces. Black circles denote triangles, and links join triangles that have non-empty intersections. Rounded-rectangles group triangles belonging to the same bloc. The right part of figure 4 shows these blocs, and the links between intersecting blocs. Each bloc has two triangulations: triangles before the vertex removal (drawn solid) are called child triangles of the bloc, and triangles after the vertex removal (in dashed) are called parent triangles.

Local analysis and synthesis

The hierarchy construction described in the previous subsection is based on the retriangulation of star-shaped polygons on the surface mesh. In this subsection, we apply the results of section 2, locally to these star-shaped polygons. This leads to a local analysis and synthesis scheme of the data set.

Let p and q be respectively the number of triangles before and after the vertex removal. Euler formula implies p = q + 2 . Let (T f i ) i=1 p denote the triangles before the vertex removal, and (T c i ) i=1 q the triangles after the vertex removal (the superscript f stands for "fine" and c for "coarse"). Since we are dealing with piecewise constant data sets, let V f denotes the space of piecewise constant functions on the triangles T f i and V c on the triangles T c i , then:

dim(V f ) = q < dim(V c ) = p = q + 2 :
The basis function ' f i in V f equals 1 on the triangle T f i and 0 elsewhere (analogously for the basis functions in V c ):

' f i = T f i ' c i = T c i :
The Gram-schmidt matrices G f and G c in V f and V c are the following diagonal matrices:

G f = diag(area(T f i ) 2 ) G c = diag(area(T c i ) 2 ) :
Therefore, equation ( 5) leads to the following local analysis matrix A:

A = 0 @ area T c i T T f j area(T c i ) 1 A
The computation of the analysis matrix B is done by first fixing a sub-matrix of B to be the identity, second solving the homogeneous linear system (6), and third performing a Gram-schmidt orthonormalization process on the rows the solution, with respect to the scalar product of matrix (G f ) ;1 . The Gram-schmidt orthonormalization process en- sures that the orthonormality condition (7) is fullfilled by the matrix B.

Global analysis and synthesis

The global analysis of a piecewise constant data set defined on an irregular planar or spherical triangulation results from the application of the local analysis formulas (1) and ( 2), with the analysis matrices described in subsection 3.2, in each bloc of the hierarchical Delaunay structure introduced in subsection 3.1, starting from the finest level, and going up to the coarsest level. After the whole analysis, each removed vertex has been assigned two wavelet coefficients. Partial reconstructions of the data set, using only wavelet coefficients of high magnitude can then be performed. This can be done by selecting a sub-hierarchy of the hierarchical Delaunay structure that contains these wavelet coefficients, and applying the reconstruction formula (4) in each bloc of this sub-hierarchy, starting from the finest level, and going up to the coarsest level. Of course, the insertion of all wavelet coefficients yields exactly the original data set. Figure 6 shows an example on a spherical data set with 1.3M faces. The data set is defined on a regular triangulation (4-to-1 split with 8 levels of subdivision starting on a icosahedron). The data value on each face derives from the ETOPO51 data set (which consists of 2160 x 4320 samples on a uniform grid). Figure 6 shows a partial reconstruction using the highest wavelet coefficients, with 150000 triangular faces (out of 1.3M). The selection of the highest wavelet coefficients clearly results in the insertion of more vertices where the data set has sharp variations, and less vertices elsewhere, thereby validating our choice for the analysis and synthesis matrices.

Figure 4 .

 4 Figure 4. Hierarchical Delaunay triangulation

Figure 5 .

 5 Figure 5. illustrates the application of the analysis scheme on one example: the top part shows the input function to the left (in the finer space V f ), the intermediate approximation in the middle (in the intermediate space Ṽ c ), and the coarser approximation to the right (in the coarser space V c ). The bottom part of figure 5 shows the two wavelet components (in the detail space W c ).

Figure 5 .

 5 Figure 5. Local analysis and synthesis: finer, intermediate and coarser approximations on top; wavelet components on bottom.

Figure 6 .

 6 Figure 6. Partial reconstructions with 150000 faces (out of 1.3M)

Available via anonymousftp at ftp://ftp.ngdc.noaa.gov/Solid Earth/Topog