\chapter{Background}
Many applications, such as Geographical Information Systems (GIS), weather
simulation and flight simulation, use terrain models in their
presentations~\cite{Li-1992-Rendering}.
Although the terrain models serves in a supporting role only, a
realistic terrain background allows the audience to visualize the results
intuitively.

Each of these applications has
different sensitivity for the accuracy of the rendered terrain models.  For
example, in GIS applications, the rendered terrain serves as the metric for
visualizing additional information, such as population density, lava flow, etc.
Hence, high accuracy terrain rendering is needed at the area of interest, but
the peripheral area can be grossly approximated.  For weather simulation, the
terrain backdrop provides the viewer with proper orientation. For example,
distinctive terrain features allow the user to determine
the direction of an in-bound cold front. The need for simulated terrain in
flight simulation is apparent.  The accuracy of the terrain used by a flight
simulation can range from tetrahedron mountains for computer gaming to
the use of highly accurate models for military flight training.

Our goal is to provide the user with real-time interactive rendering of
terrain data that is accessed remotely using a client-server system.
This thesis documents our extensive research on the methods and
algorithms used to achieve this goal.  We will first take a
closer look at the data used to model the terrain.

\section{Terrain Data}
\label{terrain-data}

Undoubtedly, the characteristics of the terrain model, such as its mathematical
properties, surface features, and digital representation, have significant
influence on rendering methods and algorithms.   A terrain data set records the
surface features of a planetary spheroid, which has a bounded surface area.
An observer, who is located outside of the spheroid, can only directly see
at most half of the surface area set at any give moment.  In addition, the 
observable surface area reduces as the observer approaching to the spheroid.
Therefore, the problem size is bounded and it is in $O(r^{2})$, where $r$ is the
radius of the spheroid.  Furthermore,
the surface of a spheroid is simple, which means an observer can go to
any point on the spheroid by walking on its surface.   This implies a unique
ordering based on the relative distance and angle from the observer.  As a
result, the rendering process can exploit this unique ordering to reduce the
amount of terrain data that it needs to render.

The characteristics of the terrain data can be used for dynamic retrieval.
The terrain features on the planetary surface can be classified by their
sizes.   The continents and oceans are the large scale terrain features.
The great plains, deserts and mountain ranges are smaller.  All
these terrain features, even the localized
terrain features, such as glaciers and lakes, are shaped by the seismic,
wind and water erosion forces.  The plate tectonics that cause a continent
to exist also introduce the pressure to form mountain ranges.  The release
of this pressure cause the earthquakes, which can divert rivers to form lakes.
As a result, the existence of similarities between various terrain features is
no surprise.  For example, one of the similarities of the terrain features are
the water/land combinations.  On the global scale, there are continents and
oceans.  When we zoom into a continent, seas and lakes are surrounded by land.
One step down in the scale, there are islands in the seas and lakes.  These
similarities offer opportunities to optimize the storage of the terrain data,
so that the detail of the terrain features can be selected at runtime.

All these terrain features are captured in the Digital Elevation Model
(DEM)~\cite{DEM}, which describe the topography of the terrain, and some
information, which is geographically related, such as satellite images and Land
Use Land Cover (LULC) data, which describes the type of vegetation.  The DEM is
a 2D array of uniformly sampled height
values of the planetary surface.  Each value represents the measurement of
height in meters relative to the sea level at that specific position.  The
satellite image are a collection of 2D array of uniformly sampled emission
values from the planetary surface.  The measurement is conducted at the
specified infrared and visible frequencies.  By connecting the points in the
DEM, a polygon based representation of the terrain is generated.  The 2D surface
images can be draped over the polygon-based DEM to create a 3D terrain model.
A standard rendering algorithm for polygon-based computer graphics objects is
then used to display the terrain model on the computer.

\section{Introduction to Computer Graphics}
\label{Intro_CG}

A polygon is formed by a collection of 3D points (vertices) and the
interconnection of these vertices (edges).  A set of polygons is used to form a
computer graphics object.  A collection of these objects forms a virtual world.
There are a wide range of rendering techniques that can be used
to generate an image of the virtual world.  Some of the techniques, such as ray
tracing~\cite{Enzmann-Ray-tracing} \cite{Wilt-Ray-tracing} and
radiosity~\cite{radiosity1}, are extremely computationally expensive, but they
can be used to
generate photo-realistic quality images.  Other polygon-based rendering
techniques, such as scan converting  Phong or Gouraud shaded polygons, are
comparatively less expensive and further speed-up can be achieved by using
graphics acceleration hardware.  Hence, these rendering methods are ideal for
interactive terrain rendering.  However, the image generated can sometimes
appear artificial to the viewer.  This is due to the approximations made by
these techniques, but we accept the trade-offs of these approximation for the
rendering speed.

\begin{figure}
\begin{center}
\begin{picture}(124,558)
\psfig{figure=standard_pipeline.eps,width=124pt,height=558pt}
\end{picture}
\end{center}
\caption{A standard eight stage rendering pipeline.  \label{standard-pipeline}}
\end{figure}

The rendering of a virtual world made up of objects represented as Phong or
Gouraud shaded polygons can be accomplished by using
a standard eight stage rendering pipeline~\cite{Foley}, which is shown in
figure~\ref{standard-pipeline}.
First, all the potentially visible objects in the virtual world
are selected.  By deselecting some visible objects, selecting some new object
from the database, and changing the orientation of
the visible objects, the virtual world appears to have a time component.

The second step is to transform all the selected objects to their proper
orientation in the virtual world.  The transformation process uses a collection
of $4 \times 4$ matrices, for translation

\begin{eqnarray}
T(x, y, z) & = &  \left[
\begin{array}{cccc}
1 & 0 & 0 & x \\
0 & 1 & 0 & y \\
0 & 0 & 1 & z \\
0 & 0 & 0 & 1 \\
\end{array} \right],
\end{eqnarray}

rotation on X, Y and Z axes

\begin{eqnarray}
R_{x}(\theta) & = &  \left[
\begin{array}{cccc}
1 & 0 & 0 & 0 \\
0 & \cos \theta & -\sin \theta & 0 \\
0 & \sin \theta & \cos \theta & 0 \\
0 & 0 & 0 & 1 \\
\end{array} \right],
\end{eqnarray}

\begin{eqnarray}
R_{y}(\theta) & = &  \left[
\begin{array}{cccc}
\cos \theta & 0 & \sin \theta & 0 \\
0 & 1 & 0 & 0 \\
-\sin \theta & 0 & \cos \theta & 0 \\
0 & 0 & 0 & 1 \\
\end{array} \right],
\end{eqnarray}

\begin{eqnarray}
R_{z}(\theta) & = &  \left[
\begin{array}{cccc}
\cos \theta & -\sin \theta & 0 & 0 \\
\sin \theta & \cos \theta & 0 & 0 \\
0 & 0 & 1 & 0 \\
0 & 0 & 0 & 1 \\
\end{array} \right],
\end{eqnarray}

and scaling (scaling is used to enlarge or reduce the size of the graphics
objects)
\begin{eqnarray}
S(s_{x}, s_{y}, s_{z}) & = &  \left[
\begin{array}{cccc}
s_{x} & 0 & 0 & 0 \\
0 & s_{y} & 0 & 0 \\
0 & 0 & s_{z} & 0 \\
0 & 0 & 0 & 1 \\
\end{array} \right].
\end{eqnarray}

%Although rotation and scaling of a 3D object can be done by using $3 \times 3$
%matrices, only $4 \times 4$ matrices can incorporate the translation operation.
The advantage of using $4 \times 4$ matrices is to combine all rotation and
scaling with translation operations into one by multiplying
all these $4 \times 4$ matrices together to form one single matrix
\begin{eqnarray}
U & = & T(x, y, z) R_{x}(\theta) R_{y}(\theta) R_{z}(\theta)
   S(s_{x}, s_{y}, s_{z}).
\end{eqnarray}
The combined matrix, $U$, is applied to each vertex of all the objects.  As
a result, the transformation of one vertex requires 16 floating-point
multiplication and 12 floating-point addition. 

Third, the objects are checked to see if they will appear in the image.
The image is the only window for the viewer to peek into the virtual world.
Similar to a real window, the image can only show part of the virtual
world.  Therefore, some of the objects in the virtual world may not appear in
the image.  Eliminating these unseen objects reduces the total amount of
rendering effort.  An object is visible if its volume intersects or is contained
within the visible volume of the virtual world.  The volume intersection
calculations can be CPU intensive, if the object has a complex shape.  A
larger transparent bounding volume with a simple shape, such as a box or
a sphere, can be used to replace the complex shape for simplifying the
calculation.  However, some of the unseen objects will be marked as visible
because of the less accurate bounding volume.
The use of tighter bounding volumes provides a more accurate visible objects
list.  A trade-off has to be made between the amount of intersection
calculation work and the amount of rendering effort saved.
The goal is to minimize the overall work of rendering the image.

Fourth, the color of the objects are computed by using the properties of the
light and the surface properties of the object.  The properties of a light
source consist of its illuminating direction, position, and color.
The illuminating direction of a light source is determined by its type, such
as directional or point source.  In order to determine if a surface is
illuminated or not by a light, the direction of the surface normal and the
illuminating direction of the light source are compared by computing the
dot product of the vectors.  If the two vectors are pointing opposite each
other, then their dot product value is positive and the surface is illuminated.
However, this does not
take into account the shadows cast by the nearby objects.  Omitting shadow
is one of the many approximations to minimize rendering resources and time.
The relative position of the light and the object determines
the illumination intensity of the object.  For example, if a surface is far
away and/or not in the illuminating direction of a light source, then the
light intensity from the light source  will be less.
If the surface is illuminated, then the color of the surface is a combination
of the object's color and the color of the light, for example a red
object under a blue light appears black.

Fifth,  the objects are transformed according to the viewing specification.
For terrain, a perspective view is typical.  The perspective transformation is
\begin{eqnarray}
M(d) & = &  \left[
\begin{array}{cccc}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & \frac{1}{1 + z_{min}} & \frac{-z_{min}}{1 + z_{min}} \\
0 & 0 & -1 & 0 \\
\end{array} \right],
\end{eqnarray}
where $z_{min} \neq -1$ and $z_{min}$ is the near plane of the viewing
volume.  The perspective transformation scales the objects by a factor inversely
proportional to the distance between the viewpoint and the objects.  The visual
effect is that all vertical parallel lines converge at a distance point and
horizontal parallel lines remain parallel. 
Hence, the rendered image conforms to our perception in the real world, where
the distant objects appear smaller but not deformed.

Sixth, the polygon based objects are clipped according to the view volume.  Any
invisible objects which escaped from step three will be eliminated entirely.
In addition, objects that are behind, too close or too far from the viewer are
also eliminated.  Because of the numerical properties of the transformation,
the objects behind the viewer become mirror reflected.   In addition,
removing the objects that are too close will avoid numerical
problems in the next stage of the pipeline.  Furthermore, objects that are far
away do not noticeably contribute to the
image, so removing them will reduce the amount of rendering effort.
Partially visible objects will have new borders that flush to the side of the
view volume.  New polygons are formed
by removing the invisible portion of the original polygon.  The
clipping computation is done by using line and plane intersection calculations.

Seventh, incorporate the viewing requirements into the objects' coordinates.
After the view transformation,
the vertex coordinates of an object changes from $(x, y, z, 1)$ to the form of
$(x', y', z', w)$.  Hence, the view transformation process can be considered as
mapping the vertices in subspace 1 onto other 3D subspaces in the 4D universe.
The $w$ coordinate of the vertices is being used to keep track of their
relative distance from the viewpoint.
This information is incorporated back to the object's 3D coordinates by dividing
the $x$, $y$ and $z$ values with the $w$ values.  Then the 3D objects are
projected onto the 2D screen space in the preparation of their rasterization.

Finally, the 2D objects are drawn onto the image buffer pixel by pixel, which
is called rasterization.  The rasterization of the 3D objects is the same as
taking a photograph of the objects, which is a parallel projection of the
objects onto a 2D space.  Then the 2D space is scaled to the desired
resolution.  However, the 2D
objects may overlap each other on the 2D plane, for example a distant 3D
object could be partially hidden behind a closer 3D object.  An arbitrary order
of rasterizing the 2D objects may put the distant objects in front of a closer
object.  The correct order of rasterization can be retrieved by sorting the list
of visible objects based on the relative distance from the viewpoint, i.e. the
$z$ coordinates of the objects.   On the other hand, a
Z-buffer records the depth information along with the color value for each 
pixel.  Only the pixel closest to the projection plane is stored, and the pixel 
contributed by objects further away will be ignored.  As a result, a correct 
image can be generated even with an arbitrary rasterization order.

The design of the standard pipeline makes no assumption on the
characteristics of the rendering objects. In order to speed up the rendering,
a partial implementation usually adapts the pipeline to exploit the objects'
characteristics. The Virtual Reality Modeling Language
(VRML)~\cite{pesce-1995-VRML} browsers, such as Webspace 1.1 VRML browser from
Silicon Graphics Inc (SGI)~\cite{Webspace} and VRweb from the Institute for
Information Processing and Computer Supported New Media (IICM), Graz University
of Technology, Austria~\cite{VRweb}, are believed to have the truest
implementation of the standard pipeline because the VRML objects are
polygon-based objects.  Both
of these browsers use OpenGL~\cite{OpenGL} for the rendering.  OpenGL provides
the basic functions, such as transformation, lighting calculation and
rasterization.  The use of OpenGL allows a portable interface to any
rendering acceleration hardware, for example SGI's RealityEngineII Graphics
pipeline and the graphics acceleration card for personal computers from
Symmetric~\cite{Symmetric} using the 3DLab GLINT graphics
processor~\cite{GLINT}.  However, specialized hardware is not the only way to
speed up computer graphics rendering.  In recent years, there have been
attempts to exploit the parallelism in computer graphics.

\section{Parallelism in Computer Graphics}

Parallel computers provide a very cost effective ways to solve
highly computational intensive problem with built-in
parallelism~\cite{gcf-Parallel-Computing} \cite{gcf-Solving-Concurrent}.
There are four types of parallelism that can be exploited: time, large grain
(task), medium grain (sub-task and/or small pocket data) and fine grain (data).
Time parallelism exists when a program is run multiple times with
different parameters.  For example, in a flight simulation program
successive images are rendered with a slight change of viewing locations.
Although these images are related in time, their rendering processes are
independent of each other, and can therefore be executed in parallel.

The different granularity of parallelism is
classified by the amount of data and the complexity of calculation
involved per parallel task.  Let us assume that the flight simulation virtual
world allows other aircraft.  Depending on the complexity of the aircraft
model and the desired detail, it may be composed of hundreds to millions of
polygons.  Since the aircraft are independent of each other, they can be
rendered separately onto separate intermediate image buffers.   Since the
amount of data and work is considerable,  this can be considered as large grain
parallelism.  The rendering of each polygon can be considered as the medium
grain parallelism. Each polygon is a self-contained data package, so all the
polygons within the virtual world can be rendered separately. 
When a polygon is rasterized, each pixel that is covered by the polygon
is visited.  However, the computation incurred by each pixel is
independent of each other.  Consequently, the work can be performed
in parallel.  Although the amount of work involve is small, the number
of pixels is large, so the exploitation of the fine-grain parallelism is
significant.

Two examples of parallel rendering that exploit the aforementioned
parallelism are a terrain rendering program from the Jet Propulsion Laboratory
(JPL), National Aeronautics and Space Administration (NASA)~\cite{Li-1996-RIVA}
and a general parallel rendering algorithm from the NASA Institute for Computer
Applications in Science and Engineering (ICASE)~\cite{Crockett-1991-rendering}.
The JPL program uses ray casting and exploits time,
task and data parallelism to obtain speedups on its rendering of terrain
data.  The program is designed to run on Multiple Instruction Multiple Data
(MIMD) parallel computers.  The $N$ nodes of the parallel computer are divided
into $M$ partitions, where $M \leq N$.  Each partition renders an image using
a full dataset with slightly different viewing parameters.  Hence, the
time parallelism is achieved.  Within each partition, the dataset is
distributed among the nodes.  Each node casts one ray per terrain datum
toward the viewer.  If the ray impacts on the image buffer then the distance
and color is compared with the buffered value.  The closest terrain data is
retained.  Otherwise, the data point is not visible.  Since each ray cast is
independent from the others, the rendering effort for the data on each node can
be considered as task parallelism.  The image buffer is also distributed among
the nodes to exploit the data parallelism associated with the pixel.

The parallel rendering algorithm from ICASE uses polygon rendering and
capitalizes on the independence of polygon and pixel.  All the polygons in the
virtual world are evenly distributed onto a MIMD computer's nodes.  All the
nodes process their polygons by using a shortened (only step two to seven)
standard rendering pipeline.  Then the polygons are redistributed according to
their screen coordinates to the node that manages that part of the image buffer.
Finally, the nodes concurrently and independently rasterize their collection of
the polygons.  The final image is formed by combining the distributed image
buffer.

\section{Multicomputers}

By distributing the work among many computing nodes, the rendering process
can be speeded up.  However, a partial image is formed on each node
after the concurrent rendering.  In order to obtain the final image,
all the intermediate image buffers must be combined in the correct order.
This is accomplished by using the communication network that connects
all the computation nodes. 

\subsection{The Message Passing Interface standard}
\label{MPI}

The communication networks of parallel computers are different.  For
example, a cluster of Alpha workstations can be 
connected by a FDDI network~\cite{NPAC}, and the software interface to
a communication channel is the standard UNIX workstation communication
library.  On the SP2 from IBM, the processing nodes are connected via a
high-performance crossbar and Ethernet.  A Message Passing Library (MPL)
provides the software interface to both networks for communicating among the
nodes. Consequently, to implement a portable program for the two parallel
computers requires a common communication library.

Message Passing Interface (MPI)~\cite{MPI} is specifically designed to
answer this need.
To ensure portability between different parallel computers, MPI provides
a logical representation of the parallel computer and a full featured
communication library, which operates on the logical machine.  Once the
logical machine and the communication library is ported to a new platform, all
parallel programs using MPI exclusively can be executed without
modification after recompilation. The porting of the logical machine is
optimally done by the manufacturer of the parallel machine. This relieves the
user from wasting any effort to customize his/her program for a specific
machine.  Hence, the user can concentrate on making algorithmic improvements to
the program.

The concern of portability of MPI program's is irrelevant, if there are no
parallel MPI program.  To encourage the user adoption of the MPI standard, MPI
provides a flexible programming environment by using an ideal architecture as
logical machine and a powerful communication library.  The logical machine
consists of a set of heterogeneous computational nodes, and they are fully
connected by a communication network. This ideal architecture allows the
programmer to easily combine algorithms, which are specifically designed for a
mesh or hypercube connected parallel computer, to form new parallel
programs.  Since these algorithms assume a specific underlying connectivity
of the computing nodes, the communication library must be capable of emulating a
given organization of nodes.  This a accomplish by using communicators.  A
communicator forms an independent logical machine with a subset of nodes.
Hence, the message sent using a communicator can not be received by other
communicators. A mesh connected parallel computer can be emulated by using a
number of row and column communicators.  Because of this flexibility and
portability, we have decide to implement our parallel rendering program using
MPI.

\section{Motivation and Project History}

The proliferation of the World Wide Web allows computer users to have instant
access to large quantities of information.  This information is typically
encoded in
HyperText Markup Language (HTML)~\cite{HTML} and is presented exactly as
in a book.  The only advantage of reading this information on the Web instead
of a book is that the reader can access other related information by
a single mouse click.  However, some of this information, such as weather data,
population density, etc,
is correlated by its geographic relationship, but this relationship can not be
easily and clearly expressed in written language.  Furthermore, the written
description only offers a single interpretation of the data.  On the other hand,
by allowing the audience to manipulate the data interactively, the data can be
presented to the reader with high fidelity.  The purpose of this project is to
provide
tools for laying out the geographically related information for educational
application.  Students are selected as the initial audience because training can
be provide through their teachers and funding can be obtained to add new
equipment to the school, such as a new computer network and some personal
computers.  With this equipment, a client-server computing environment
can be set-up to deliver the parallel computer's computing power to the school.
Our pilot software, which consists of a client program and a number of server
side programs, render a satellite image of the Albany, New York area on top of
a 3D terrain model.

The data sets, the satellite image and the DEMs, are transformed by our
preprocessing program into a 2D grid of polygon based graphics objects.
A polygon based representation of the data is chosen because computer
hardware and software is readily available for rendering polygon based computer
graphics objects.  Since this is the popular techniques used for engineering
rendering and some computer games,
Gouraud shaded images of these objects are rendered using
a custom software parallel rendering engine on an 8 node IBM SP2 parallel
computer at Northeast Parallel Architectures, Syracuse University. The high
speed communication network of the SP2 is used by the
rendering engine to assemble the resulting images.  To reduce the rendering
time of these images,  a subset of the terrain objects are selected according
to the visibility dictated by the viewing parameters.  After the images are
completed, they are shipped across the network to our client program, which
display the images as they arrive. The client program also provides the user
control to adjust the viewing parameters.

With the accelerated demand of computer rendering capability fueled by computer
games and the
use of the Web, the mass production of specialized graphics hardware has been
set in motion.  As a result, the use of a parallel computer for rendering
graphics for the general public is no longer cost effective, although it would
still be useful for applications with very large data sets that require very
high quality rendering and only have a single user at a time.  Even though
development of the custom programs was suspended before they were fully
completed, they give us an in-depth look into what rendering techniques are
needed for rendering terrain.  This proved to be valuable
knowledge for the design and implementation of Virtual Reality Modeling
Language (VRML)~\cite{pesce-1995-VRML} terrain rendering system.

Our VRML terrain rendering system uses a number of VRML built-in
capabilities, such as texture map and level of detail (LOD), to provide an
interactive presentation of the terrain data.  Because VRML browsers must
process the entire VRML file before rendering the first image, the amount of
memory and CPU power of client machine are the limiting factors of the
complexity of the VRML world and the frame rate of the rendering. The VRML
system uses texture mapped polygons, so that the number
polygons in the terrain graphics objects can be minimized.  However, the
terrain graphics objects of the original terrain data may
contain too many polygons for some of the clients. A real-time wavelet based
terrain model optimization algorithm was designed to reduce the number polygons
in the graphical objects to a more manageable number, while retaining the
dominant features of the real estate.  To increase the quality of the
presentation, the client's rendering power must be efficiently used.  Since
the VRML browser is not optimized for terrain rendering, the designer of the
virtual world must give hints to the browser on when and where it can conserve
the rendering effort.  This can be achieved by providing multiple levels of
detail for the same object and instructions on when to display what level of
detail of the object.  By using our optimization algorithm a multi-resolution
terrain model can easily generated. As a result, the browsers can provide high
quality images with better frame rate by minimizing the number of visible
polygons in the terrain model based on the viewing parameters.

\section{Overview of this Thesis Contents}

In chapter~\ref{viewing-characteristics}, we explore the viewing
characteristics of the
terrain data.  By exploiting these viewing characteristics, we can reduce
amount of the terrain data and minimize the rendering
effort.  A real-time wavelet-based terrain data reduction algorithm is
presented in chapter~\ref{multi-terrain}.  In
chapter~\ref{parallel-terrain-rendering}, we present our design of a
parallel terrain rendering algorithm.  A web-based terrain rendering
system using VRML is described in chapter~\ref{VRML_terrain}.  In
chapter~\ref{terrain-rndering-systems}, the server-side and client-side
rendering approaches are compared and we describe future research
directions for terrain rendering.


