From alvin@npac.syr.edu Mon Sep 16 20:43:39 1996
Date: Thu, 12 Sep 1996 14:35:50 -0400 (EDT)
From: Alvin Leung <alvin@npac.syr.edu>
To: Paul Coddington <paulc@beetlejuice.npac.syr.edu>
Subject: New Chapter 3


\chapter{Multiresolution Terrain Model}

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

The justification of using a simplified model is based on the assumptions that
less information requires less processing.
A computer graphical 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 have to be multiplied by at least one 4x4 transformation
matrix.  When the object is finally rendered on screen, the connectivity
information provided by the polygon is used to draw the object.  As a result,
an object
with smaller number of vertices and polygons (low resolution model) usually
renders faster than the same object with larger number of vertices and polygons
(high resolution model).  Chapter~\ref{rendering} provides more detail
on computer graphics rendering.

The lower resolution model is constructed by sacrificing 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 the
low resolution model.  A naive method of generating low resolution model by
uniformly subsampling the original model can ``miss'' the Grand Canyon, if none
of the sampling points lie within the canon.  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 turning parameters for controlling
the amount of detail used in constructing the model.

We briefly present a number of existing mesh optimization algorithms.
In~\cite{Schroeder-1992-decimation}, the article introduces a general purpose
mesh reduction algorithm.  It eliminates vertices that are within a user
specified tolerance from the replacement polygon.  It also retains sharp
transition 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 object with curve
surfaces.  It can be applied to terrain because at the current resolution of
the terrain data, sharp transitions are
seldom because the sparse sample of the height values only capable of
capturing large and slow varying terrain features.  However, there exists
large terrain features and detail terrain model, where curve surfaces is not
a good fit for the model.  In the paper ``Mesh
Optimization''~\cite{hoppe-1993-mesh}, a brute force optimization techniques
are used to find a better placement for the new set of vertices.  The
algorithm generates very good simplified model, but the runtime is in the
order of minutes for a small set of input points.  Scape~\cite{scape-paper}
implements an algorithm specifically designed for optimizing terrain grid.  It
has a short runtime and generates accurate simplified terrain mesh.  Some of
these algorithms~\cite{turk-1992-retiling}~\cite{hoppe-1993-mesh} is
intended for generating static simplified model.  All of them assumes the user
have access to the full data set and powerful computer to handle the
optimization program.

In this chapter, we will describe a wavelet based algorithm for constructing
a multiresolution terrain model.  This algorithm not only generates a lower
resolution terrain model with the dominant features, but it is also able to
construct the model on-the-fly with an adjustable parameter for
the amount of detail.  This ``dial-a-terrain-model'' feature is
used to generate VRML terrain, and the details of its use will be describe in
chapter~\ref{VRML_terrain}.  In the next section, we will cover the
mathematical background of this algorithm.

\section{Introduction to Wavelet}

In this section, we will briefly introduce the wavelet transformation, and
these information is obtained from~\cite{1994-SIGGRAPH-CourseNotes}
and~\cite{1995-SIGGRAPH-CourseNotes}.  Wavelet transformation is derived
from 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(t)$ is the window function.
If the window size is selected to vary according to the spread of frequency
component 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 have relatively long smooth area between localized
discontinuities.  As our description of the terrain data on
section~\ref{terrain-data},  we believe that the terrain data have similar
characteristics.  Hence, transformation similar to Windowed Fourier
transformation is suitable for analyze and model our terrain data.

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, Wavelet transformation uses 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 normization 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 not chosen at the
ideal size, the time and frequency resolution for the scaled window is lower
than proper selected window size.
According to~\cite{1994-SIGGRAPH-CourseNotes}, the advantage of using Wavelet
transformation is that the user uses any bandpass function for $h(t)$.  A
bandpass function satisfy three properties,
\begin{itemize}
\item regular
\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 discreteize version of the transformation is used to work with the numerical
data.

Using Wavelet transformation to model terrain requires to apply the
transformation to many different scale (see section~\ref{terrain-data} for the
different scale 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 part of
the output, the input signal of the filter can be recovered.  Hence, the
different scale of terrain can be peeled off layer by layer, and recombine as
needed.  We will discuses in more detail on
wavelet transformation and the filter bank 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 function, 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 generated without loss of accuracy.
However, the polygonal representation of the sparse (DEM) may introduce some
inaccuracy as presented in section~\ref{terrain_sparse_error_discuession}.

\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.
\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 operates on a cluster of four
adjacent height values, and the box functions do not overlap each other.
The box function can trivially apply 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-arraies with power of 2 dimension.  The drawback of
the non-overlapping box function is that any height variation among the
clusters can not be captured, and rapidly varying of height is possible
at the boundary of the 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 can not faithfully capture any terrain features that is less than
200m in diameter (see
section~\ref{terrain-data}), so the terrain cross section form the DEM is
usually very
smooth.  As a result, the discontinuity between the box function do not
necessary induce any loss of information because there is no information to
lose.  On the other hand, this lack of tie 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 inner product 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 of $V^{i}$ and
$W^{i}$.  Using a normalized basis allows us to conduct a meaningful comparison
on the numerical coefficient obtained from different scales.  This is an
important requirement for the 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 put a upper bound of on the amount of storage requirement.
The non-overlapping condition indicates that the information encoded by the
wavelet coefficients are only from one box.  This limits the complexity of the
computation because the recovery computation can not 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 on 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 on 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, standard decomposition scheme
(figure~\ref{standard_decomposition})
and non-standard decomposition scheme as shown in
figure~\ref{non-standard_decomposition}.  We elected to use the non-standard
decomposition for two reasons.  First, the results of each iteration of the
decomposition has 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 dimension.  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 are
\[
\begin{array}{ll}
 & \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] \\
\end{array}.
\]
The $\frac{1}{\sqrt{2}}$ is the normalization factor for the bases.
For larger terrain model, the transformation matrices becomes
\[
\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].
\]
Rewrite the matrix form into equations
\[
\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}.
\]
A couple of remarks arise from the above equations.  First,
$T^{i - 1}(x + 0.5, y + 0.5)$ is two times of 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 difference 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 raggedness 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
magnetite of the coefficients not only compares the coefficients' effect on the
reconstructed height values, but also takes into account of the area affect by
the coefficients.

The computation for reconstructing the terrain model from the
wavelet coefficients are given by
\[
\begin{array}{ll}
 & \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}.
\]
The computation is expressed in equations form,
\[
\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 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)$ is equal to zero, then
\[
T^{i}(x, y) = T^{i}(x + 1, y) = T^{i}(x, y + 1) = T^{i}(x + 1, y + 1).
\]
The coplanar characteristic of these points and the polygonal representation of
the terrain model enables us to store a sparse representation of the terrain
model 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 wavelet terrain is divided into the decomposition and
reconstruction stage.  The decomposition phase converts the terrain model into a
collection of the wavelet coefficients.  Then a subset of the wavelet
coefficients is chosen and used to reconstruct a terrain model.  Since only
input for the reconstruction phase is the wavelet coefficients, the
decomposition need only be done once.

   As shown in the previous section, the wavelet terrain model can be
constructed by using sparse matrix multiplication.  Although the sparse matrix
multiplication offers advantages (as stated in
~\cite{Lounsbery-1994-Topological}), such as, code modularity, ease of
programming
by using existing software library, and ease of changing the basis functions,
we choose not to use it.  Form our past
experience~\cite{LU-1990}, a fast sparse matrix operation requires a careful
exploitation of the sparsity pattern, such as, 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 algorithms is designed
based on the use of this code.

\subsection{Decomposition of 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 the following:

\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}

The number of operations is given by a geometric series with its 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 the all 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 use a database to help us implement this accessing scheme because it
provides us a flexible the data storage management.  As a result, we do not
need
to design, build, test and maintain a set of subroutines to manage this
elaborate storage scheme.  Furthermore, we can flexibly select the coefficient
by using a single SQL command.  The grid is implemented as a collection of
tables.  The name of the table identifies a grid element.  Each entry of the
tables consists of the followings:

\begin{tabbing}
123\=coefficients\=1\=123\=123\=123\=123\=123\=123\= \kill
\>k \> \>level identification \\
\>x, y \> \>the coordinates of the box within the level \\
\>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 according to
the descending order of k and is stored into a link list.  The reconstruction
stage uses this collection of coefficients to form the terrain model.

\subsection{Reconstruction of 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.  Then,
a triangulation algorithm is used to generate the connectivity of the polygons.

The reconstruction algorithm transforms $M$ coefficient entries, where
$M \leq \lfloor \frac{N^{2}}{3} \rfloor + 1$, in the link list, $L$, which is
sorted in descending order of k, into a terrain model.  Since the
decomposition algorithm recusively devided the size of DEM array by 2,
the total number of wavlet entries is equal to the geometric series, where the
constant equal to $\frac{N^{2}}{4}$ and the ratio equal to $\frac{1}{4}$, plus
the last elevation entry.  The algorithm is initialized with the last
elevation value, the first entry of $L$.  While the link 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 exists, 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.
The following is the above algorithm in pesudo code for the reconstruction of 
a terrain model from the link list, $L$:

\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}

After the reconstruction of the height\_values(k, i, j), it needs to be
mapped
into the terrain model's coordinates system.  There are two possible mapping.
First, the height\_values(k, i, j), where $N = 2^{k}$, is from the original
terrain grid.  In this
case, the i and j is multiplied by the separation factor between two data points
to put the vertex in its proper location.
Otherwise, the height\_values(k, i, j), where $2^{k} < N$, is from a later
layer.  It is multiplied by $\frac{2^{k}}{N}$ to
compensate for the factor of 2 difference between two adjacent layers.
The i and j is needed to be enlarged by a factor of 2 for each layer.  As a
result, the i and j is multiplied by $\frac{N - 1}{2^{k} - 1}$ and the
separation factor.

\begin{figure}
\begin{center}
\begin{picture}(369,387)
\psfig{figure=terrain-vertices.eps,width=369pt,height=387pt}
\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}(369,387)
\psfig{figure=terrain-polygon.eps,width=369pt,height=387pt}
\end{picture}
\end{center}
\caption{After triangulation is completed, the sparse DEMs become terrain
model.  Each line represents an edge of a polygon.
\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 model of 90Km NW of Albany. 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 model of 90Km NW of Albany. 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
have dramatic impacts on the complexity of the reconstruction algorithm.
Assuming the retrieval of height values is $O(1)$, the computation
complexity of the reconstruction algorithm is $O(M)$.
Our implementation stores the height values in a binary search tree.
Initially, we assume that
the seemingly ``random'' insertion order and repeated deletions makes a
reasonable balanced binary search tree.  The runtime of this implementation
is shown in table~\ref{terrain-runtime-nobalance}.  A second implementation
uses AVL binary tree, and it runtime is at least 4 times better than the
previous version (table~\ref{terrain-runtime-balance}).  As a result, the
access time of the data play a significant
role on 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 model.  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,
the vertices is projected onto the X-Y plane and a triangulation is performed
on the 2D points.  This process enable us to form polygons by connecting the
vertices.  There are two criteria that the triangulation must compiles.
\begin{enumerate}
\item The vertices must be connected to its nearest neighbors.
\item The triangle can not overlap.
\end{enumerate}
According to~\cite{lee-1980-two}, the circumcircle of any triangle in the
Delaunay triangulation contains no other vertex in its interior.  In addition,
a typical
implementation of the algorithm avoids overlapping triangles.  Consequently,
the two criteria is fulfilled.  We selected to use
a stable implementation of the Delaunay triangulation algorithm,
triangle~\cite{triangle}, from Jonathan Shewchuk at Carnegie Mellon University.
The result is shown in figure~\ref{terrain-polygon}.
Since Delaunay triangulation algorithm is a well
studied area~\cite{guibas-1985-primitives}, we quote from the
reference~\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 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 apply 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.  We show the output terrain models generated
by this algorithm.  In the next section, we will presents our results to
support our claim in the previous sections.

\section{Findings}
We have given a detailed account of the wavelet terrain algorithm previously.
Now, we present some results to support our claim of
automatic generation of sparse terrain models and 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,229)
\psfig{figure=terrain_sparse_error.eps,width=432pt,height=229pt}
\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 (left) and the spares (right) 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 coincide with a
transition between two DEM data points as shown in
figure~\ref{terrain_sparse_error}.  These transitions usually occur at
the boundary of the water bodies, which is represented by a cluster equal
height DEM data points.  The magnitude of the error depends on the
differences of the two points.  In our empirical studies, these errors are not
detectable in 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 model
originally contains $16,384$ vertices and $32,258$ triangles.  The Crater lake
model originally contains $226 \times 459$ vertices and requires $306,860$
polygons to model.
Since our algorithm can only handle terrain model 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 lefthand 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 exploit 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 on number of vertices and polygons \label{terrain-reduction}}
\end{table}

\label{terrain-multi-resolution}
Our algorithm can also be used to generate a high quality simplified terrain
model.  In figures~\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 with 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} listed the runtime of the generation of
various sparse models.  As a result, we have accomplished designing an
algorithm to generate 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 table shows the number of wavelet coefficients, vertices and
polygons of different model 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 table shows the timing of different stages for generating
different model 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 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,228)
\psfig{figure=terrain_accuracy.eps,width=432pt,height=228pt}
\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 have only one adjustable parameter, which is the choice of
the wavelet coefficients used in reconstruction phase.
Hence, the selection policy of the wavelet coefficients have significant
impact on the quality of the terrain model.  Our current selection is made via
a user adjustable threshold.  Any wavelet coefficient entry, which have 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.  The higher values implies that the least 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 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 give a good result for
terrain visualization, it
has two disadvantages. First, the total number of vertices and polygons in the
model can not 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 user's familiarity with the terrain contour.  Keeping a
statistical model
for each data block may be able to give hint for the user.  Second, the cut
off value can not be used to compute the numerical accuracy of the terrain
model.  As shown is figure~\ref{terrain-accuracy}, the tolerance factor
determines the magnitude of error between the new vertex relative to the
original vertices.  However, additional error can be introduced by the
placement of
the new vertices, when the edges join the new vertex with the neighboring
vertices.  As a result, the current coefficient selection is only suitable for
visualizing terrain data.



Alvin - a devil with shades
  ------------------------------------------------------------------------------
mailto:alvin@nova.npac.syr.edu
www: http://www.npac.syr.edu/users/alvin


