\chapter{Multiresolution Terrain Model}
\label{multi-terrain}

A multiresolution terrain model consists of a collection of simplified, but
``similar looking'' sets of polygons created using Digital Elevation Models
(DEMs) and terrain image data.  As stated in
section~\ref{high2low-flight}, the viewing requirements and characteristics
for terrain rendering
dictate that a mix of different resolution terrain
models is essential to reduce the required computing power for a satisfactory
viewing experience.

The justification for using a simplified model is based on the assumption that
less information requires less processing.
A computer graphics object is defined by a collection of 3D points (vertices)
and the connectivity of the points (polygons) in a 3D virtual world.
To place the object at a specific location within the virtual world, 
each of the vertices has to be multiplied by at least one $4 \times 4$
transformation matrix.  When the object is finally rendered on screen, the
connectivity information provided by the polygons is used to draw the object.
As a result, an object
with a smaller number of vertices and polygons (low resolution model) usually
renders faster than the same object with a larger number of vertices and
polygons (high resolution model).  Section~\ref{Intro_CG} provides more detail
on computer graphics rendering.

The lower resolution model is constructed by sacrificing some of the modeling
detail.  However, if the wrong information is discarded, then the topography of
the terrain
can be changed drastically.  Imagine if the Grand Canyon was omitted from a
low resolution model for Arizona.  A naive method of generating a low resolution
model by
uniformly subsampling the original model can ``miss'' the Grand Canyon, if none
of the sampling points lie within the canyon.  A ``similar looking'' low
resolution model must retain the dominant features of the terrain.  Hence, a
more sophisticated algorithm is needed to produce the low resolution model.
Furthermore, this algorithm must have some tuning parameters for controlling
the amount of detail used in constructing the model.

First, we briefly present a number of existing mesh optimization algorithms,
which might be used for creating multiresolution terrain models.
Schroeder, Zarge and Lorensen~\cite{Schroeder-1992-decimation} provide a
general purpose
mesh reduction algorithm.  It eliminates vertices that are within a
user-specified tolerance from the replacement polygon.  It also retains sharp
transitions of the model.  The algorithm essentially acts as a filter that
removes redundant information and noise from the model.  Truk's
algorithm~\cite{turk-1992-retiling} is designed to handle 3D objects with curved
surfaces.  It can be applied to terrain because at the current resolution of
the terrain data, sharp transitions are
rare because the sparse sampling of the height values is only capable of
capturing large and slowly varying terrain features.  However, there exist
large terrain features and detailed terrain models, where curved surfaces are
not a good fit for the model.  Hoppe {\sl et al}.~\cite{hoppe-1993-mesh}
present a mesh optimization strategy using general purpose optimization
techniques, such as steepest descent or simulated annealing.  The goal of the
optimization method is to find a better placement for the new set of vertices.
The algorithm generates a very good simplified model, but the runtime is on the
order of minutes for a small set of input points.  Scape~\cite{scape-paper}
implements an algorithm specifically designed for optimizing terrain grids.  It
has a short runtime and generates an accurate simplified terrain mesh.  Some of
these algorithms~\cite{turk-1992-retiling}~\cite{hoppe-1993-mesh} are
intended for generating a simplified persistent data model, which has a
fixed complexity.  If different resolution is desired, then the original data
has to be processed again.  Therefore, all of these algorithms assume the user
has access to the full data set and a powerful computer to handle the
optimization program once for each level of complexity needed.

In this chapter, we will describe a wavelet-based algorithm for constructing
a multiresolution terrain model.  This algorithm not only generates lower
resolution terrain models that retain the dominant features, but is also
able to construct the model on demand with an adjustable parameter for
the amount of detail.
This ``dial-a-terrain-model'' is made possible by storing the model's 
wavelet coefficients, and the graphical object is reconstructed when needed.
This feature is ideal for serving graphical objects on the web because the
complexity of the model can be adjusted to suite the client's needs and
capabilities.  We use this feature to generate VRML terrain, and the details of
its use will be described in
chapter~\ref{VRML_terrain}.  In the next section, we will cover the
mathematical background of this algorithm.

\section{Introduction to Wavelets}

In this section, we briefly introduce the wavelet
transformation~\cite{1994-SIGGRAPH-CourseNotes}
\cite{1995-SIGGRAPH-CourseNotes} \cite{rioul-wavelet}.
The wavelet transformation is derived
from the Windowed Fourier transformation with scaled window size.  The Windowed
Fourier transformation is formulated as
\begin{equation}
   F_{w}(\tau, f) = \int_{-\infty}^{\infty} f(t)w^{\ast}(t - \tau)
   e^{-2i\pi ft}\,dt, \label{Windowed-Fourier-transform}
\end{equation} where $\tau$ is the position of the filter, $f(t)$ is the
signal and $w^{\ast}(t)$ is the convoluted window function.
If the window size is selected to vary according to the spread of frequency
components of the signal, this allows higher resolution in time at the
discontinuities, and
retains a good frequency resolution at high frequency.  The scheme only
works when the signal has a relatively long smooth area between localized
discontinuities.  As our description of the terrain data in
section~\ref{terrain-data} implies, we believe that terrain data has similar
characteristics.  Hence, a transformation similar to the Windowed Fourier
transformation is suitable to analyze and model our terrain data.

The wavelet transformation is derived by defining 
the basic wavelet as $h(t) = w(t)e^{-2i\pi ft}$.  To simplify the selection of
window size, the wavelet transformation uses a scaled window size, and
the scaled basic wavelet becomes $h_{a} = \frac{1}{\sqrt{|a|}}h(\frac{t}{a})$,
where $a$ is the scale and $\frac{1}{\sqrt{|a|}}$ is the normalization factor.
The wavelet transformation is formulated as
\begin{equation}
   WF(\tau, a, f) = \frac{1}{\sqrt{|a|}}\int_{-\infty}^{\infty}
   f(t)h^{\ast}(t - \tau)\,dt. \label{Wavelet-transform-forward}
\end{equation}
The scaled window size provides similar properties as selecting the window size
according to the frequency spread.  Since the window size is a fixed scale, it
cannot be chosen at the ideal size for every function, the time and frequency
resolution for the scaled window is lower than the optimal resolution.
According to~\cite{1994-SIGGRAPH-CourseNotes}\cite{rioul-wavelet}, the advantage
of using the wavelet transformation is that the user can use any bandpass
function for $h(t)$, as long as it satisfies the following three properties:
\begin{itemize}
\item regular~\cite{rioul-wavelet}
\item $\int h^{2}(t)\,dt$ is finite
\item $\int h(t)\,dt = 0$
\end{itemize}
The inverse Wavelet transformation is given by
\begin{equation}
   f(t) = \int_{-\infty}^{\infty} WF(\tau, a, f)
      \frac{1}{\sqrt{|a|}}h(\frac{t - \tau}{a})\,\frac{da\,dt}{a^{2}}.
      \label{Wavelet-transform-back}
\end{equation}
A discretized version of the transformation is used to work with the numerical
data.

Using the wavelet transformation to model terrain requires applying the
transformation to many different scales (see section~\ref{terrain-data} for the
different scales of terrain).  A filter bank can be setup by using a series of
low-pass and high-pass filters.  After a signal passes through one set of low
and high-pass filters, the signal is divided into two parts.  The output of
the low-pass filter is the original signal minus some details and it continues
its journey down the filter bank.  By combining the two corresponding parts of
the output, the input signal of the filter can be recovered.  Hence, the
different scales of terrain can be peeled off layer by layer, and recombined as
needed.  We will discuss the
wavelet transformation and the filter bank in more detail in the next section.

\section{Mathematical Background}
\label{wavelet-box}
\label{sparse-terrain}

The computation of the multiresolution terrain model is based upon the
definition of a set of scaling functions, an inner product defined on the
scaling functions and the wavelet for the orthogonal complement space.
Lounsbery's thesis~\cite{Lounsbery-1994-Topological} and SIGGRAPH 94 Course
Notes~\cite{1994-SIGGRAPH-CourseNotes} provide details on the
mathematical derivation.  The scaling functions span a vector space and the
function also partitions the vector space into a nested vector space,
\[
V^{0} \subset V^{1} \subset \cdots \subset V^{j - 1} \subset V^{j},
\]
where $j$ denotes the relative ranking of the vector space, based on
the different scale.  This mathematical feature allows us to extract
similarity across different scales of the terrain data.  As a result, a sparse
Digital Elevation Model (DEM) can be losslessly generated by eliminating
duplicated parts with the same height, such as lakes, rivers and plains.
However, the polygonal representation of the sparse DEM may introduce some
inaccuracy, as explained in section~\ref{terrain_sparse_error_discuession}.
Furthermore, we can use the lossless sparse DEM to generate a good lossy
representation as our simplified (lower resolution) model.

\begin{figure}
\begin{center}
\begin{picture}(257,331)
\psfig{figure=box.eps,width=257pt,height=331pt}
\end{picture}
\end{center}
\caption{Two different scales of box functions operate on a 2D array.
The missing connections are the result of using a non-overlapping box function.
\label{box_function}}
\end{figure}

The box functions are chosen as the scaling functions because they
are a convenient choice for operating on a 2D array of height values, as shown
in figure~\ref{box_function}.  Each box function averages on a cluster of four
adjacent height values, and the box functions do not overlap each other.
The box function can trivially be applied to a 2D array with power-of-2
dimension.
For an array with other than power-of-2 dimension, the array can be divided
into a collection of sub-arrays with power-of-2 dimension.  The drawback of
the non-overlapping box function is that any height variation among the
clusters cannot be captured, and rapid variation of height is possible
at the boundary of two adjacent boxes.  During the construction of the low
resolution DEM, some terrain characteristics may be over attenuated because of
the missing information.
However, this limitation does not pose any significant impact on our terrain
model because of the resolution of the current DEM.
The 100m DEM cannot faithfully capture any terrain features that are less than
200m in diameter (see section~\ref{terrain-data}), so the terrain cross section
from the DEM is usually very
smooth.  As a result, the discontinuity between the box functions does not
necessarily induce any loss of information because there is no information to
lose.  On the other hand, this lack of any ties between the boxes allows the
independent recovery of each individual box.  As a result, a subsection of the
terrain model can be recovered with work proportional to the sum of all boxes
enclosed within the
subsection at each level.  A detailed description of recovering a subsection of
the terrain is available in section~\ref{terrain-model-subsection}.

The existence of an inner product in the vector space allows us to find all the
orthogonal functions in $V^{i}$ of
the functions in $V^{i - 1}$, where $0 < i \leq j$.  These
orthogonal functions form a orthogonal complement space, $W^{i - 1}$. The
orthogonal complement space allows us to capture the loss information for the
process of mapping $V^{i}$ to $V^{i - 1}$.  Hence, all the functions in
$V^{i}$ can be recovered without loss.  In addition, the inner product also
provides a normalization factor for both basis sets of $V^{i}$ and
$W^{i}$.  Using a normalized basis allows us to conduct a meaningful comparison
on the numerical coefficients obtained from different scales.  This is an
important requirement for runtime tunability.  The standard inner product
is used to construct a orthogonal complement space, $W^{j}$.

The Haar wavelets are used as the basis of $W^{i}$ because they
are orthogonal to the box function. They also provide the two essential
conditions, compact support and being non-overlapping, to make the recovery of
each box a constant amount of work.  A wavelet is said to have compact support,
when its computation involves a finite number (in practice, a small number) of
coefficients.  This puts an upper bound on the amount of storage required.
The non-overlapping condition indicates that the information encoded by the
wavelet coefficients is only from one box.  This limits the complexity of the
computation because the recovery computation cannot be proportional to the
number of all wavelet coefficients.  As a result, the recovery of each box
function in $V^{i}$ involves a small number of coefficients for the box
functions in $V^{i -1}$ and for the Haar wavelets in $W^{i - 1}$.  In the next
section, we will discuss the numerical considerations for preparing a wavelet
terrain model.

\section{Numerical Aspects}

\begin{figure}
\begin{center}
\begin{picture}(400,305)
\psfig{figure=standard_decomposition.eps,width=400pt,height=305pt}
\end{picture}
\end{center}
\caption{Standard decomposition of a 2D array.
$T^{i}$ and $C^{i}_{0}$ represents a matrix.
\label{standard_decomposition}}
\end{figure}

\begin{figure}
\begin{center}
\begin{picture}(400,396)
\psfig{figure=non-standard_decomposition.eps,width=400pt,height=396pt}
\end{picture}
\end{center}
\caption{Non-standard decomposition of a 2D array. $T^{i}, C^{i}_{0},$
$C^{i}_{1}$ and $C^{i}_{2}$, where $0 < i \leq j$, represents a matrix.
\label{non-standard_decomposition}}
\end{figure}

After a brief discussion of the components of the wavelet terrain model, we now
describe the various numerical aspects of the decomposition and
reconstruction of the terrain model.  There are two ways to apply the box
function to a 2D array of height values, a standard decomposition scheme
(shown in figure~\ref{standard_decomposition})
and a non-standard decomposition scheme (shown in
figure~\ref{non-standard_decomposition})~\cite{1994-SIGGRAPH-CourseNotes}
\cite{1995-SIGGRAPH-CourseNotes}, where $T^{a}$ is a 2D array of height
values, $C^{a}_{b}$ is a 2D array of wavelet coefficients, $a$ is the
iteration index of the decomposition and $b$ is the unique identifier of the
wavelet coefficient array.  We elected to use the non-standard
decomposition for two reasons.  First, the results of each iteration of the
decomposition have the expected shapes and sizes that are implicitly suggested
by the concept of the nested vector space.  Therefore, we can verify the
correctness of the
implementation by rendering the height values in $T^{i}$, where $0 < i \leq j$.
Second, we believe that the non-standard decomposition scheme allows the wavelet
computation to extract similarity across different scales of the terrain data
in 2 dimensions.  As a result, a better sparse representation of the terrain can
be generated.

The numerical computation in matrix form for decomposing the DEM into wavelet
coefficients is
\[
\begin{array}{lcl}
\left[
\begin{array}{cc}
T^{i - 1}& C^{i - 1}_{0}\\
C^{i - 1}_{2}& C^{i}_{1}\\
\end{array}
\right]
& = & 
cQT^{i}P \\
\end{array},
\]
where $c$ is a scalar normalization factor, $Q$ and $P$ is the transformation
matrices for the bases.
The transformation for the simplest 2 $\times$ 2 terrain model is
\begin{equation}
\begin{array}{lcl}
cQT^{i}P
& = & 
\frac{1}{\sqrt{2}}
\left[
\begin{array}{cc}
1 & 1 \\
1 & -1 \\
\end{array}
\right]^{T}
\left[
\begin{array}{cc}
T^{i}(x, y) & T^{i}(x, y + 1) \\
T^{i}(x + 1, y) & T^{i}(x + 1, y + 1) \\
\end{array}
\right]
\frac{1}{\sqrt{2}}
\left[
\begin{array}{cc}
1 & 1 \\
1 & -1 \\
\end{array}
\right] \\
& = & 
\frac{1}{2}
\left[
\begin{array}{cc}
1 & 1 \\
1 & -1 \\
\end{array}
\right]^{T}
\left[
\begin{array}{cc}
T^{i}(x, y) & T^{i}(x, y + 1) \\
T^{i}(x + 1, y) & T^{i}(x + 1, y + 1) \\
\end{array}
\right]
\left[
\begin{array}{cc}
1 & 1 \\
1 & -1 \\
\end{array}
\right] \\
\end{array}.
\label{2by2_decomposition}
\end{equation}
The transformation matrices, $Q$ and $P$, are the a transpose of each other for
the Haar wavelet, and the $\frac{1}{2}$ is the normalization factor for the
bases (see reference~\cite{Lounsbery-1994-Topological} for more detail).
For the general case, the transformation matrix becomes
\begin{equation}
\begin{array}{lcl}
P & = &
\left[
\begin{array}{llllllllll}
1 & 0 & \cdots & 0 & 0 &  1 &  0 & \cdots & 0 & 0 \\
1 & 0 & \cdots & 0 & 0 & -1 &  0 & \cdots & 0 & 0 \\
0 & 1 & \cdots & 0 & 0 &  0 &  1 & \cdots & 0 & 0 \\
0 & 1 & \cdots & 0 & 0 &  0 & -1 & \cdots & 0 & 0 \\
0 & 0 & \cdots & 0 & 0 &  0 &  0 & \cdots & 0 & 0 \\
\vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots
& \vdots \\
0 & 0 & \cdots & 0 & 0 &  0 &  0 & \cdots & 0 & 0 \\
0 & 0 & \cdots & 1 & 0 &  0 &  0 & \cdots &  1 &  0 \\
0 & 0 & \cdots & 1 & 0 &  0 &  0 & \cdots & -1 &  0 \\
0 & 0 & \cdots & 0 & 1 &  0 &  0 & \cdots &  0 &  1 \\
0 & 0 & \cdots & 0 & 1 &  0 &  0 & \cdots &  0 & -1 \\
\end{array}
\right] .
\end{array}
\end{equation}
Rewriting the matrix form of equation~\ref{2by2_decomposition} into equations
gives
\[
\begin{array}{clc}
T^{i - 1}(x + 0.5, y + 0.5) & = & \frac{T^{i}(x, y) + T^{i}(x + 1, y)
+ T^{i}(x, y + 1) + T^{i}(x + 1, y + 1)}{2} \\
\\
C^{i - 1}_{0}(x + 0.5, y + 0.5) & = & \frac{T^{i}(x, y) + T^{i}(x + 1, y)
- T^{i}(x, y + 1) + T^{i}(x + 1, y + 1)}{2} \\
\\
C^{i - 1}_{1}(x + 0.5, y + 0.5) & = & \frac{T^{i}(x, y) - T^{i}(x + 1, y)
- T^{i}(x, y + 1) - T^{i}(x + 1, y + 1)}{2} \\
\\
C^{i - 1}_{2}(x + 0.5, y + 0.5) & = & \frac{T^{i}(x, y) - T^{i}(x + 1, y)
+ T^{i}(x, y + 1) - T^{i}(x + 1, y + 1)}{2} \\
\end{array}.
\]
The $(x, y)$ coordinate indicates the placement of the terrain elevation entry,
but the coordinate for the wavelet coefficients is there for uniformity.
A couple of remarks arise from the above equations.  First,
$T^{i - 1}(x + 0.5, y + 0.5)$ is two times the average height of the four
points enclosed by the box function.  Second,
$C^{i - 1}_{0}(x + 0.5, y + 0.5), C^{i - 1}_{1}(x + 0.5, y + 0.5)$ and
$C^{i - 1}_{2}(x + 0.5, y + 0.5)$ are the differences between the height value,
$T^{i - 1}(x + 0.5, y + 0.5),$ at level $i - 1$ and its corresponding height
values, $T^{i}(x, y), T^{i}(x + 1, y), T^{i}(x, y + 1),$ and
$T^{i}(x + 1, y + 1)$, within the box at level $i$.
As a result, the magnitude of the coefficients is a good indicator of the
terrain's jaggedness within the box.
Furthermore, they are also scaled up by a factor of 2.
The scaling factor of 2 between each layer reflects the extended influence
of the height value and the wavelet coefficients, which cover
a bigger piece of real estate in the later layer.  Consequently, a selection
based on the magnitude of the coefficients not only compares the effect of the
coefficients on the reconstructed height values, but also takes into account
the area affected by the coefficients.

The computation for reconstructing the terrain model from the
wavelet coefficients is given by
\begin{equation}
\begin{array}{lcl}
T^{i} & = & \frac{1}{\sqrt{2}}
\left[
\begin{array}{cc}
1 & 1 \\
1 & -1 \\
\end{array}
\right]
\left[
\begin{array}{cc}
T^{i - 1}(x + 0.5, y + 0.5) & C^{i - 1}_{0}(x + 0.5, y + 0.5) \\
C^{i - 1}_{2}(x + 0.5, y + 0.5) & C^{i - 1}_{1}(x + 0.5, y + 0.5) \\
\end{array}
\right]
\frac{1}{\sqrt{2}}
\left[
\begin{array}{cc}
1 & 1 \\
1 & -1 \\
\end{array}
\right]^{T} \\
\end{array}.
\label{2by2_reconstruction}
\end{equation}
The computation of equation~\ref{2by2_reconstruction} expressed in equation
form is:
\[
\begin{array}{lcl}
T^{i}(x, y) & = &
\frac{T^{i - 1}(x + 0.5, y + 0.5) + C^{i - 1}_{2}(x + 0.5, y + 0.5)
+ C^{i - 1}_{0}(x + 0.5, y + 0.5) + C^{i - 1}_{1}(x + 0.5, y + 0.5)}{2} \\
\\
T^{i}(x, y + 1) & = &
\frac{T^{i - 1}(x + 0.5, y + 0.5) + C^{i - 1}_{2}(x + 0.5, y + 0.5)
- C^{i - 1}_{0}(x + 0.5, y + 0.5) + C^{i - 1}_{1}(x + 0.5, y + 0.5)}{2} \\
\\
T^{i}(x + 1, y) & = &
\frac{T^{i - 1}(x + 0.5, y + 0.5) - C^{i - 1}_{2}(x + 0.5, y + 0.5)
+ C^{i - 1}_{0}(x + 0.5, y + 0.5) - C^{i - 1}_{1}(x + 0.5, y + 0.5)}{2} \\
\\
T^{i}(x + 1, y + 1) & = &
\frac{T^{i - 1}(x + 0.5, y + 0.5) - C^{i - 1}_{2}(x + 0.5, y + 0.5)
- C^{i - 1}_{0}(x + 0.5, y + 0.5) - C^{i - 1}_{1}(x + 0.5, y + 0.5)}{2} \\
\end{array}.
\]
Each of the above equations can be computed by only 4 floating point additions
and one multiplication.  The simplicity of the computation makes this algorithm
suitable for implementation on a less powerful computer (see
chapter~\ref{VRML_terrain} for more detail).
Finally, we would like to point out a special case. If 
$C^{i - 1}_{0}(x + 0.5, y + 0.5), C^{i - 1}_{1}(x + 0.5, y + 0.5)$ and
$C^{i - 1}_{2}(x + 0.5, y + 0.5)$ are all equal to zero, then
\begin{equation}
T^{i}(x, y) = T^{i}(x + 1, y) = T^{i}(x, y + 1) = T^{i}(x + 1, y + 1).
\label{DEM_sparse}
\end{equation}
As a result, we can exploit this characteristic and maintains a
sparse representation of the DEM without loss of accuracy.  We will describe
the algorithm and implementation of the wavelet terrain model in the next
section.

\section{Algorithm and Implementation}
\label{low-overhead-access}

   The implementation of a wavelet terrain representation is divided into a
decomposition and a
reconstruction stage.  The decomposition phase converts the terrain model into a
collection of wavelet coefficients.  Then a subset of the wavelet
coefficients is chosen and used to reconstruct a terrain model.  Since the only
input for the reconstruction phase is the wavelet coefficients, the
decomposition need only be done once in a preprocessing step.

As shown in the previous section, the wavelet terrain model can be
constructed by using sparse matrix multiplication.  Although sparse matrix
multiplication offers advantages (as stated in
~\cite{Lounsbery-1994-Topological}) such as code modularity, ease of
programming by using existing software libraries, and ease of changing the
basis functions, we chose not to use it.  From our past
experience~\cite{LU-1990}, a fast sparse matrix operation requires a careful
exploitation of the sparsity pattern, that is, the distribution and density of
the non-zero values.  The high overhead of a general sparse matrix operation
subroutine on memory and accessing data limits its typical performance to the
collection of its design targeted matrices.  During the terrain model
reconstruction phase, the
wavelet coefficients are chosen dynamically, so that the sparsity pattern of
the wavelet coefficient matrix varies dramatically.  As a result, we decided
to implement a low overhead data accessing code and the algorithm is designed
based on the use of this code.

\begin{figure}
\begin{tabbing}
123\=123\=123\=123\=123\=123\=123\=123\= \kill
\>for (k = j; k $>$ 0; k$++$) \\
\>\>for (x = 0; x $< 2^{k - 1}$; x$++$) \\
\>\>\>for (y = 0; y $< 2^{k - 1}$; y$--$) \\
\>\>\>\> Compute $T^{k - 1}(x, y)$, $C^{k - 1}_{0}(x, y)$,
$C^{k - 1}_{1}(x, y)$ and $C^{k - 1}_{2}(x, y)$ \\
\>\>\>\> Store $T^{k - 1}(x, y)$ \\
\>\>\>\> if $C^{k - 1}_{0}(x, y)$ or $C^{k - 1}_{1}(x, y)$
or $C^{k - 1}_{2}(x, y) > 0$ then \\
\>\>\>\>\> Store $C^{k - 1}_{0}(x, y)$,
$C^{k - 1}_{1}(x, y)$ and $C^{k - 1}_{2}(x, y)$ \\
\end{tabbing}
\caption{The decomposition algorithm in pseudo code form.
\label{decompsition_alg}}
\end{figure}

\subsection{Wavelet Decomposition of the Digital Elevation Model}
The decomposition algorithm transforms a $N \times N$ array of height values,
where $N = 2^{j}$ and $j = 1, 2, 3, \cdots$, to a collection of wavelet
coefficients.  The algorithm is shown in figure~\ref{decompsition_alg}.
The number of operations is given by a geometric series with a constant equal
to $N^{2}$ and the ratio between each adjacent term equal to $\frac{1}{2^{2}}$.
Hence, the computation complexity is $O(N^{2})$.  The memory accessing pattern
is shown in figure~\ref{decomposition-memory}, and the memory requirement is
also $O(N^{2})$.

\begin{figure}
\begin{center}
\begin{picture}(312,72)
\psfig{figure=decomposition-memory.eps,width=312pt,height=72pt}
\end{picture}
\end{center}
\caption{Memory access pattern for decomposition of wavelet terrain
\label{decomposition-memory}}
\end{figure}

The low overhead data accessing scheme (previously mentioned in
section~\ref{low-overhead-access}) uses two indexing methods to achieve fast
access on the coefficients.  The first indexing is accomplished by imposing a 2D
grid onto the entire terrain model, and the decomposition is performed within
each of the grid elements.  This grid forms a 2D hash table for all wavelet
coefficients.  The second indexing is done by processing all the wavelet
coefficients within the grid element.  As a result, the
work of accessing a subset of wavelet coefficients consists of a single
transformation from the geographical coordinates to the grid coordinates plus
the work needed to locate the data entries within the grid element.
Since all numerical computations are done with the coefficients inside
a box (referring to the box described in section~\ref{wavelet-box}), we can
further reduce the number of data entries by grouping all related
coefficients together.  Consequently, only one high overhead access per box
is needed during reconstruction.

We avoided the need to design, build, test and maintain a set of subroutines to
manage this elaborate storage scheme by
using a database to help us implement this accessing scheme, because it
provides us with flexible data storage management.  We can select a set of
coefficients from the database by using SQL commands.  The grid
is implemented as a collection of tables in the database.  The name of the
table identifies a grid element.  Each entry of the tables consists of the
following:

\begin{tabbing}
123\={\tt coefficients}\=1\=123\=123\=123\=123\=123\=123\= \kill
\>{\tt k} \> \>level identification \\
\>{\tt x, y} \> \>the coordinates of the box within the level \\
\>{\tt coefficients} \> \>$C^{k}_{0}(x, y), C^{k}_{1}(x, y)$ and
$C^{k}_{2}(x, y)$
\end{tabbing}

After a query is performed, the result is ordered by the database in descending
order of {\tt k} and is stored into a linked list.  The reconstruction stage
uses the collection of coefficients to form the terrain model.

\subsection{Reconstruction of the Terrain Model}
\label{reconstruction}

The reconstruction of a terrain model is divided into two steps.  First, the
height values are constructed from the wavelet coefficients.  Next,
a triangulation algorithm is used to generate the connectivity of the polygons.

The reconstruction algorithm transforms $M$ coefficient entries
in the linked list, $L$, which is
sorted in descending order of {\tt k}, into a terrain model.  Since the
decomposition algorithm recursively divided the size of the DEM array by 2,
the total number of wavelet entries is equal to the geometric series,
$a\sum_{i = 1}^{j}r_{i}$, where $a$ is $\frac{N^{2}}{4}$ and $r$ is
$\frac{1}{4}$, plus the elevation array $T^{0}$, which has the size of 1, so
that $M \leq \lfloor \frac{N^{2}}{3} \rfloor + 1$.  The algorithm is initialized
with the last
elevation value, the first entry of $L$.  While the linked list is not empty,
the next set of wavelet coefficients, $C_{0}$, $C_{1}$, $C_{2}$ and the
corresponding elevation, $T^{i}(x, y)$, is used to compute the new height
values, $T^{i + 1}(x, y)$, $T^{i + 1}(x, y + 1)$, $T^{i + 1}(x + 1, y)$ and
$T^{i + 1}(x + 1, y + 1)$.  Since the algorithm is also used to reconstruct
a sparse DEM, $T^{i}(x, y)$ may not exist because of the DEM sparse
representation, so a substitute value for
$T^{i}(x, y)$ is used.  The substitution value is selected from the set of
$T^{j}(x, y)$, where $j = 0, \ldots, i - 1$, and the value of $j$ is largest.
Figure~\ref{reconstruction_alg} presents the above algorithm in pseudo code for
the reconstruction of a terrain model from the linked list, $L$.

\begin{figure}
\begin{tabbing}
123\=123\=123\=123\=123\=123\=123\=123\= \kill
\>[T, k, i, j] = list\_head(L) \\
\>delete\_list\_head(L) \\
\>height\_values(k + 1, i, j) = T \\
\>while (!is\_empty(L)) \\
\>\>[$C_{0}$, $C_{1}$, $C_{2}$, k, i, j] = list\_head(L) \\
\>\>delete\_list\_head(L) \\
\>\>if height\_values(k, i, j) exists then \\
\>\>\>T = height\_values(k, i, j) \\
\>\>\>delete height\_values(k, i, j) \\
\>\>else \\
\>\>\>k' = k - 1 \\
\>\>\>i' = i / 2 \\
\>\>\>j' = j / 2 \\
\>\>\>while (height\_values(k', i', j') not exists) then \\
\>\>\>\>k' = k' - 1 \\
\>\>\>\>i' = i' / 2 \\
\>\>\>\>j' = j' / 2 \\
\>\>\>T = height\_values(k', i', j') \\
\>\>\>if height\_values(k', i', j') occupies the same location in the terrain
model as \\
\>\>\>\>height\_values(k + 1, 2i, 2j), height\_values(k + 1, 2i + 1, 2j), \\
\>\>\>\>height\_values(k + 1, 2i, 2j + 1) and
height\_values(k + 1, 2i + 1, 2j + 1) then \\
\>\>\>\> Tag height\_values(k', i', j') \\
\>\>compute height\_values(k + 1, 2i, 2j), height\_values(k + 1, 2i + 1, 2j), \\
\>\>\>height\_values(k + 1, 2i, 2j + 1) and
height\_values(k + 1, 2i + 1, 2j + 1) \\
\>\>\>with T, $C_{0}$, $C_{1}$, $C_{2}$ \\
\> delete tagged height\_values \\
\end{tabbing}
\caption{The reconstruction algorithm in pseudo code form.
\label{reconstruction_alg}}
\end{figure}

After the reconstruction of the array {\tt height\_values(k, i, j)}, it needs
to be mapped
into the terrain model's coordinate system.  There are two possible mappings.
We recall that $N = 2^{l}$, where $l = 1, 2, 3, \cdots$, is the dimension of
the 2D array of height values.
If $N = 2^{k}$, then {\tt height\_values(k, i, j)} is from the original
terrain grid.  In this
case, {\tt i} and {\tt j} are multiplied by the separation factor between
two data points to put the vertex in its proper location.
Otherwise, if $2^{k} < N$, {\tt height\_values(k, i, j)} is from a later
layer, in which case it is multiplied by a factor of
$\frac{2^{k}}{N} = 2^{k - l}$.
This compensates for the normalization factor of 2, the difference between two
adjacent layers, which is introduced at the decomposition time.
We also recall that every decomposition step from $T^{i}$ to $T^{i - 1}$
reduces the array dimension by a factor of 2.  Since {\tt i} and {\tt j} are
indexing the height value in the array of $T^{i - 1}$, a scaling factor for the
separation factor of {\tt i} and {\tt j} is needed.
The scaling factor is computed by dividing the number of separations in the
original 2D height value array which is ($N - 1$) by the number of spaces
between data elements of the array for the current height value, which is
($2^{k} - 1$).  As a result, {\tt i} and {\tt j} are multiplied by
$\frac{N - 1}{2^{k} - 1}$ and the separation factor.

\begin{figure}
\begin{center}
\begin{picture}(432,474)
\psfig{figure=terrain-vertices.eps,width=432pt,height=474pt}
\end{picture}
\end{center}
\caption{This figure shows the density of the vertices for the sparse DEM with
different coefficient cutoff values.  Each black dot represents a height
value of the DEM. \label{terrain-vertices}}
\end{figure}

\begin{figure}
\begin{center}
\begin{picture}(432,474)
\psfig{figure=terrain-polygon.eps,width=432pt,height=474pt}
\end{picture}
\end{center}
\caption{After triangulation is completed, the sparse DEMs become a terrain
model.  Each line represents an edge of a polygon.  The models are generated
with progressively increasing error tolerance.  There is a collection of long
skinny triangles on right hand side of the three most sparse models.
These are the result of the forced boundaries of the sparse models.
\label{terrain-polygon}}
\end{figure}

\begin{table}
\begin{tabular}{|l|l|l|l|} \hline
model threshold & database access & wavelet transformation & triangulation \\
\hline
dense & 1.15 sec & 1.86 sec & 1.76 sec \\ \hline
0 & 1.38 sec & 1.67 sec & 1.66 sec \\ \hline
5 & 0.713 sec & 1.13 sec & 1.01 sec \\ \hline
10 & 0.570 sec & 0.830 sec & 0.670 sec \\ \hline
20 & 0.309 sec & 0.230 sec & 0.361 sec \\ \hline
40 & 0.162 sec & 0.0670 sec & 0.177 sec \\ \hline
80 & 0.0680 sec & 0.0159 sec & 0.072 sec \\ \hline
\end{tabular}
\caption{The table shows the timing of different stages for generating
different resolution terrain models for an area of $128 \times 128$ data
points. The
elevations are stored in a binary tree.  \label{terrain-runtime-nobalance}}
\end{table}

\begin{table}
\begin{tabular}{|l|l|l|l|} \hline
model threshold & database access & wavelet transformation & triangulation \\
\hline
dense & 1.12 sec & 0.409 sec & 1.69 sec \\ \hline
0 & 1.17 sec & 0.380 sec & 1.60 sec \\ \hline
5 & 0.740 sec & 0.238 sec & 0.989 sec \\ \hline
10 & 0.547 sec & 0.167 sec & 0.675 sec \\ \hline
20 & 0.333 sec & 0.0895 sec & 0.372 sec \\ \hline
40 & 0.167 sec & 0.0421 sec & 0.172 sec \\ \hline
80 & 0.0956 sec & 0.0284 sec & 0.073 sec \\ \hline
\end{tabular}
\caption{The table shows the timing of different stages for generating
different resolution terrain models of an area of $128 \times 128$ data points.
The elevations
are stored in a balanced binary tree.  \label{terrain-runtime-balance}}
\end{table}

The distribution pattern of height values for a sparse terrain
model is no longer uniform; hence, a flexible storage and retrieval mechanism
is needed.  In addition, the access time of storage and retrieval
has dramatic impacts on the complexity of the reconstruction algorithm.
Assuming the retrieval of height values is $O(1)$, the computational
complexity of the reconstruction algorithm is $O(M)$.
Our implementation stores the height values in a binary search tree.
Initially, we assumed that
the seemingly ``random'' insertion order and repeated deletions would make a
reasonably balanced binary search tree.  The runtime of this implementation
is shown in table~\ref{terrain-runtime-nobalance}.  A second implementation
(table~\ref{terrain-runtime-balance}) uses an AVL binary tree~\cite{Algorithm},
and its runtime is at least 4 times better than the previous version.  All
timings are done on a 4 processor Silicon Graphics (SGI) ONYX.  The access 
time of the data plays a significant role in slowing down the algorithm.  A
faster accessing scheme, such as a hash table with a
balanced binary search tree for bucket management, is recommended for the
generation of larger terrain models.  Our estimated runtime of
our fastest implementation is $O(M\log n)$, where $n$ is
the number of height values in the binary tree.
A set of outputs from the wavelet decompression is shown in
figure~\ref{terrain-vertices}.

After the computation of the vertices and their height values,
the vertices are projected onto the X-Y plane and a triangulation is performed
on the 2D points.  This process enables us to form polygons by connecting the
vertices.  Because the original DEM is a digital representation of a simple
surface, the sparse DEM also inherits the same properties.
Hence, each data point has an implicit connectivity to its nearest
neighbor.  Furthermore, a graphical reconstruction of the simple surface does
not involve overlapping triangles, and maximization of the triangle's angles 
reduces computer graphics artifacts, while rendering the surface.
Therefore three criteria that the triangulation must satisfy are:
\begin{enumerate}
\item Each vertex must be connected to its nearest neighbors.
\item The triangles cannot overlap.
\item All angles of the triangles should be maximized
\end{enumerate}

These criteria are satisfied by Delaunay triangulation, which is a unique
triangulation of a set, $S$, of
$n$ points in $d$ dimensional Euclidean space (for details on the
mathematical definition, proofs and properties see~\cite{fortune-1992-voronoi}
\cite{lee-1980-two}), having two important properties:
\begin{enumerate}
\item there is no element of S inside the circumsphere of any triangle.
In 3D case, a circumsphere is the sphere, which the three vertices of a triangle
lies on its surface.
\item In 2D, Delaunay triangulation maximize the minimum angle of the
resulting triangles.
\end{enumerate}
Since there are no other points inside the circumsphere of any triangle, the
vertices of the triangle are the closest elements in $S$ to each other.  Hence,
each vertex is connected to its nearest neighbors.  Furthermore, this
guarantees that there are no overlapping triangles, because overlapping
triangles implies connection to neighbor further away than necessary.  The
angles of the triangles are maximized when Delaunay triangulation is performed
in 2D.  We chose to use a stable implementation of the Delaunay triangulation
algorithm, {\tt triangle}~\cite{triangle}, from Jonathan Shewchuk at Carnegie
Mellon University.  The result is shown in figure~\ref{terrain-polygon}.
On the right edge of the three out of four triangulated DEMs in figure
~\ref{terrain-polygon}, there is a collection of long skinny triangles.  These
are the result of the forced boundaries of the sparse models.  In
chapter~\ref{VRML_terrain}, we propose a corrective measure to counter this
problem.  Since the Delaunay triangulation algorithm is a well
studied area~\cite{guibas-1985-primitives}, we merely quote from Guibas, Knuth
and Sharir~\cite{guibas-1992-randomized} that the complexity of
the Delaunay triangulation algorithm is $O(n\log n)$, where $n$ is
the number of vertices.
Since $M$ coefficient entries is less than $n$ vertices, the overall runtime
for the reconstruction algorithm is bounded by the runtime of the
Delaunay triangulation algorithm, which is $O(n\log n)$.

\label{terrain-model-subsection}
  Cropping a subsection of the terrain model is an operation that needs to be
done very often. For example, the whole terrain data set is divided into a
number of pieces before a wavelet transformation is applied to the data.  In order
to expand a
subsection of the terrain model, we need to expand all boxes that contain
the subsection.  This will produce the $T^{0}(x,y)$ for the subsection.
Then the coefficients, which are from the boxes overlapped or contained within
the subsection, are selected.  Finally, height values which are outside of the
subsection boundary, are removed.  In the next section, we show the output
terrain models generated by this algorithm and we present our results.

\section{Findings}
In this chapter, we have given a detailed account of the wavelet terrain
algorithm.
Now, we present some results to support our claim that wavelets can be used
for automatic generation of sparse terrain models and for generating high
quality multi-resolution terrain
models with a single adjustable parameter.  We also show the
capabilities and the shortcoming of our current controlling scheme for the
parameter.

\begin{figure}
\begin{center}
\begin{picture}(432,237)
\psfig{figure=terrain_sparse_error.eps,width=432pt,height=237pt}
\end{picture}
\end{center}
\caption{The inaccuracy of the sparse terrain model
\label{terrain_sparse_error}}
\end{figure}

\begin{figure}
\begin{center}
\begin{picture}(222,288)
\psfig{figure=crater_dense.ps,width=222pt,height=288pt}
\end{picture}
\
\begin{picture}(222,288)
\psfig{figure=crater_000.ps,width=222pt,height=288pt}
\end{picture}
\end{center}
\caption{The dense, 306,860 polygons (top), and the sparse, 281,473 polygons
(bottom), terrain model of Crater Lake area is shown.  \label{crater_sparse_0}}
\end{figure}

\label{terrain_sparse_error_discuession}

In section~\ref{sparse-terrain}, we claimed that our wavelet terrain algorithm
can convert a dense terrain model into a sparse terrain model without adding any
significant error.  These errors
only occur where the missing linkage between the two boxes coincides with a
transition between two DEM data points.  These transitions usually occur at
the boundaries of water bodies, which are represented by a cluster of equal
height DEM data points, so that one box (for the water) can be represented by
a single midpoint value, while a neighboring box stays at the original values
(see figure~\ref{terrain_sparse_error}).  The magnitude of the error depends on
the differences between the two points.  In our empirical studies, these errors
are not detectable in a terrain visualization application
(figure~\ref{crater_sparse_0}).
Table~\ref{terrain-reduction} shows the percentage of reduction on the number
of vertices and polygons for three sparse terrain models.  The two smaller
models originally both contain $16,384$ vertices and $32,258$ triangles.  The
Crater Lake
model, which is bundled with Scape software~\cite{scape}, originally contains
$226 \times 459$ vertices and requires $306,860$ polygons to model.
Since our algorithm can only handle terrain models with size of power-of-2, we
padded the model with zero values and the final size is $512 \times 512$
vertices.  The Crater Lake data is stored in the section starting from the
upper left-hand corner $(0, 0)$ to the lower right corner $(335, 458)$.  Even
with the enlarged model, the sparse Crater Lake model contains less vertices
and polygons than the dense Crater Lake model.  Although the amount of reduction
depends on the terrain's contour enclosed within the block, the algorithm 
automatically exploits the smoothness and generates sparse data, when possible.

\begin{table}
\begin{tabular}{|l|l|l|} \hline
format and location &  vertices (\% reduction) & polygons (\% reduction) \\
\hline
dense 90Km NW of Albany & 16384 & 32258 \\ \hline
sparse 90Km NW of Albany & 14836 (9.4\%) & 29203 (9.5\%) \\ \hline
dense Albany & 16384 & 32258 \\ \hline
sparse Albany & 10693 (34.7\%) & 21023 (34.8\%) \\ \hline
dense Crater Lake & 262144 & 522242 \\ \hline
sparse Crater Lake & 141143 (44.9\%) & 281473 (46.1\%) \\ \hline
\end{tabular}
\caption{Reduction in number of vertices and polygons for dense and sparse
terrain models. \label{terrain-reduction}}
\end{table}

\label{terrain-multi-resolution}
Our algorithm can also be used to generate a high-quality simplified (lower
resolution) terrain
model.  In figure~\ref{crater_model_compare},
we show a sequence of successively simplified terrain models with the dense
terrain model as reference.
The sparse model shows
similar terrain features as the dense version, but uses less polygons.
The number of vertices, polygons and the number of
wavelet coefficients used in the reconstruction of the sparse models
are listed in table~\ref{terrain-size-crater}.
Table~\ref{terrain-runtime-crater} lists the run times for the generation of
various sparse models.  As a result, we have accomplished our goal of designing
an algorithm to generate a multiresolution terrain model.

\begin{table}
\begin{tabular}{|l|l|l|l|} \hline
model threshold & wavelet coefficients & vertices & polygons \\
\hline
0 & 47048 & 141143 & 281473 \\ \hline
5 & 20778 & 62351  & 124238 \\ \hline
10 & 12304 & 36927 & 73483 \\ \hline
20 & 6817 & 20471 & 40673 \\ \hline
40 & 3386 & 10176 & 20157 \\ \hline
80 & 1898 & 5709 & 11282 \\ \hline
160 & 1064 & 3196 & 6290 \\ \hline
320 & 696 & 2091 & 4101 \\ \hline
640 & 507 & 1522 & 2984 \\ \hline
\end{tabular}
\caption{The number of wavelet coefficients, vertices and
polygons of different resolution models of Crater Lake
\label{terrain-size-crater}}.
\end{table}

\begin{table}
\begin{tabular}{|l|l|l|l|} \hline
model threshold & database access & wavelet transformation & triangulation \\
\hline
dense & 19.71 sec & 13.47 sec & 44.34 sec \\ \hline
0 & 10.50 sec & 3.86 sec & 24.51 sec \\ \hline
5 & 4.44 sec & 1.70 sec & 8.91 sec \\ \hline
10 & 2.82 sec & 1.64 sec & 5.15 sec \\ \hline
20 & 1.52 sec & 0.53 sec & 2.35 sec \\ \hline
40 & 0.79 sec & 0.28 sec & 0.99 sec \\ \hline
80 & 0.51 sec & 0.14 sec & 0.53 sec \\ \hline
160 & 0.29 sec & 0.10 sec & 0.30 sec \\ \hline
320 & 0.21 sec & 0.06 sec & 0.27 sec \\ \hline
640 & 0.14 sec & 0.05 sec & 0.14 sec \\ \hline
\end{tabular}
\caption{The timing of different stages for generating
different resolution models of Crater Lake. The elevations are stored in a
balanced binary tree.  \label{terrain-runtime-crater}}
\end{table}

\begin{figure}
\begin{center}
\begin{picture}(140,181)
\psfig{figure=crater_dense.ps,width=140pt,height=181pt}
\end{picture}
\
\begin{picture}(140,182)
\psfig{figure=crater_040.ps,width=140pt,height=182pt}
\end{picture}
\
\begin{picture}(140,182)
\psfig{figure=crater_080.ps,width=140pt,height=182pt}
\end{picture}
\
\begin{picture}(140,182)
\psfig{figure=crater_160.ps,width=140pt,height=182pt}
\end{picture}
\
\begin{picture}(140,182)
\psfig{figure=crater_320.ps,width=140pt,height=182pt}
\end{picture}
\
\begin{picture}(140,180)
\psfig{figure=crater_640.ps,width=140pt,height=180pt}
\end{picture}
\end{center}

\caption{The dense (upper left most) and various sparse terrain models of
the Crater Lake date set is shown.  The sparse models are generated with
threshold set to 40, 80, 160, 320 and 640, respectively.
\label{crater_model_compare}}
\end{figure}

\begin{figure}
\begin{center}
\begin{picture}(432,248)
\psfig{figure=terrain_accuracy.eps,width=432pt,height=248pt}
\end{picture}
\end{center}
\caption{The inaccuracy of the simplified terrain model
\label{terrain_accuracy}}
\end{figure}

\label{wavelet-parameter}
As indicated by our reconstruction algorithm in section~\ref{reconstruction},
wavelet terrain has only one adjustable parameter, which is the choice of
the number of wavelet coefficients used in the reconstruction phase.
Hence, the selection policy of the wavelet coefficients has significant
impact on the quality of the terrain model.  Our current selection is made via
a user-adjustable threshold.  Any wavelet coefficient entry, which has at
least one coefficient with absolute magnitude larger than the threshold, is
selected.  The cut-off value represents a detail filtering factor for the
terrain model.  A higher value implies that less detail is included in
the model.  Since the wavelet coefficients are normalized, the effect of the
threshold is halved on each successive layer.  Consequently, an effective way
to create a multiresolution model is to multiply the initial cut-off value by a
factor of a power of two.  This will ensure a uniform increase of smoothness
over the entire model.

Although this selection policy gives a good result for
terrain visualization, it
has two disadvantages. First, the total number of vertices and polygons in the
model cannot be easily predicted.  Since the wavelet coefficients' magnitude
is determined by the terrain contour, the proper selection of the threshold
by the user depends on the user's familiarity with the terrain contour.  Keeping
a history of tolerance versus the number of polygons used
for each data block may be able to give hints to the user (see
section~\ref{stat-parameter} for more detail).  Second, the cut-off
value cannot be used to compute the numerical accuracy of the terrain
model.  As shown in figure~\ref{terrain_accuracy}, the tolerance factor
determines the magnitude of error between the new vertex relative to the
original vertex.  When an edge joins a new vertex with a neighboring vertex,
the edge implicitly defines a linear interpolation function for the interval
between the two end height values.  An additional error can be introduced by
the placement of the new vertex because currently all new vertices are placed at
the center of their box function.  In the case shown in
figure~\ref{terrain_accuracy}, the
additional error is eliminated if the new vertex is placed at the same
location of the original vertex on the left side.  However, proper placement
strategy for the new vertices is computationally expensive, so we compromise
on the accuracy of the models to gain fast generation of the models.  Because
of the existence of uncontrolled error in the generated model, we do not believe
the current wavelet terrain algorithm is suitable for other purposes, such as a
ground vehicle simulator~\cite{scape-paper}, where a highly ``accurate''
simplified
model is needed.  As a result, our algorithm will have to be customized for
entertainment use, where the accuracy of the model is less important, but
the rendering speed is tightly controlled.  The client computer's
the number of polygons per second rendering
capability can be obtained or estimated dynamically.
Consequently, the number of polygons of the model is the preferred
tuning parameter.  In the next section, we will introduce a method to enable
us to select a terrain model based on the number of polygons.

\section{Improvements}
\label{stat-parameter}
  We have mentioned in the previous section that the value of the current
adjustable parameter for the wavelet terrain algorithm does not provide a
prediction for the number of polygons in the terrain model, due to the fact
that each terrain block has an unknown unique function for mapping the error
tolerance to the number of polygons used in the model.  This relationship for
the Crater Lake model is shown in figure~\ref{crater_polygon}.  This function is
highly non-linear and its characteristics cannot be easily modeled by
analytical means.  Furthermore, there are many terrain blocks and each of them
needs a prediction model, so we need to employ a low overhead and low
cost modeling method in order to model each mapping function separately.

We propose to model the mapping function by
a piecewise linear function.  The minimum initial data required by the
piecewise linear function is two data points, i.e., two reconstructions for each
terrain block. These data points can be store as a part of the
meta-data for that terrain block.  Hence, the overhead is $O(n)$, where $n$ is
the number of
terrain blocks in the entire database.  We recommend that the two data points
be selected at the beginning and the end of the typical operation 
spectrum of the error tolerance.  The estimation of the number of polygons in
the model at a given error tolerance level is given by
the following equation
\begin{equation}
\begin{array}{lcl}
y & = & \frac{y_{1} - y_{0}}{x_{1} - x_{0}}x
+ \frac{x_{1}y_{0} - x_{0}y_{1}}{x_{1} - x_{0}} \label{piecewise-linear},\\
\end{array}
\end{equation}
where $(x_{0}, y_{1})$ and $(x_{1}, y_{1})$ are two data point from the
prediction model, $y$ is the number of polygons in the terrain model and $x$ is
the error tolerance level.
Each prediction requires three floating point additions and two floating point
multiplications only, so the computational cost is very low.  In addition to low
overhead and low cost, this method also allows us to improve our
prediction accuracy over time.  Since
the prediction calculation only requires two data points; we can insert
additional data points into our prediction
model without incurring any additional computational cost during the estimation
calculation.  When there are more
than two points in the prediction model, we select the two closest data points,
preferably one before and one after the target point, for the calculation.

The advantage of this adaptive approach is that a popular terrain model will
have a good prediction assuming not all the users are requesting the model at
the same error tolerance. A new data point for the prediction model is
available after the terrain model is
reconstructed with a error tolerance value not in our prediction model.
Figure~\ref{crater-linear-predict} shows a progressive refined prediction
model for the Crater Lake model.  A table of maximum error of the
prediction model is shown in table~\ref{prediction-error}.
With two data points,  the prediction model has a discrepancy of 1860 percent.
However, as more data points are inserted, the prediction model improves.
When six data points are used for the prediction, the maximum error in various
range of error tolerance levels varies from 378 percent to 1 percent.
The prediction error between error
tolerance interval of 1 to 40 remains at 378 percent because the new terrain
models were generated using higher error tolerances, 80, 160 and 320.  With
these three additional data points, the maximum prediction error between 40 to
640 is down from 1860 percent to 12 percent.  As a
result, the prediction model provides a reliable means to provide a ``guess''
of the number of polygons used in the terrain model.  Furthermore, the
adaptiveness of the prediction model is based on the usage pattern of the
model.

\begin{figure}
\begin{center}
\begin{picture}(360,252)
\psfig{figure=crater_polygons.eps,width=360pt,height=252pt}
\end{picture}
\end{center}
\caption{The log-log plot shows the highly non-linear function between the
error tolerance and the number of polygons used in the Crater Lake terrain
model.
\label{crater_polygon}}
\end{figure}

\begin{figure}
\begin{center}
\begin{picture}(360,252)
\psfig{figure=crater_fit.eps,width=360pt,height=252pt}
\end{picture}
\end{center}
\caption{The plot shows five piecewise linear function modeling the relationship
between the error tolerance and the number of polygons for the Crater Lake
model.  A total of six terrain models is reconstructed.  The models are
generated with threshold set to 1, 40, 80, 160, 320 and 640, respectively.
\label{crater-linear-predict}}
\end{figure}

\begin{table}
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|} \hline
              & \multicolumn{5}{c|}{Error tolerance range} \\ \hline
  data points & 1 to 40 & 40 to 80 & 80 to 160 & 160 to 320 & 320 to 640 \\ \hline
  2           & 1187 & 1693 & 1860 & 1856 & 1433 \\ \hline
  3           &  378 & 43   & 67   &  71  & 59   \\ \hline
  4           &  378 & 12   & 20   &  27  & 23   \\ \hline
  5           &  378 & 12   & 4    &  8   &  8   \\ \hline
  6           &  378 & 12   & 4    &  8   &  1   \\ \hline
\end{tabular}
\end{center}
\caption{This table shows the maximum percent error between the predicted and
the actual number of polygons used in the Crater Lake model.
\label{prediction-error}}
\end{table}


