\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 given 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 terrain data is simply connected, which means no simple closed curve in the terrain data set surrounds points not in the data set.  This implies a simple coordinate system can be imposed onto the data set; for example, 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.  The seismic, wind and water erosion forces shape all terrain features, even the localized terrain features, such as glaciers and lakes.  The plate tectonics that cause a continent to exist also introduce the pressure to form mountain ranges.  The release of this pressure causes 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 is the water/land combination.  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 is 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 collection of 3D points (vertices) and the interconnection of these vertices (edges) form a polygon.  A set of polygons is used to form a computer graphics object.  A collection of these objects forms a virtual world. There are wide ranges 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-off of these approximations 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 object 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 is 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 type of light source, such as directional or point, determines the illuminating direction.  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}
\label{parallel_CG}

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 amount of data and the complexity of calculation involved per parallel task classify the different granularity of parallelism.  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 data set with slightly different viewing parameters.  Hence, the time parallelism is achieved.  Within each partition, the data set 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.  Combining the distributed image buffer forms the final image.

\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 FDDI network connects a cluster of Alpha workstations~\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 programs is irrelevant, if there are no parallel MPI programs.  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 audiences 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 rendering engine to assemble the resulting images uses the high speed communication network of the SP2.  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-rendering-systems}, the server-side and client-side rendering approaches are compared and we describe future research directions for terrain rendering.

