\chapter{Parallel Terrain rendering}
\label{parallel-terrain-rendering}

We have implemented two methods for visualizing our terrain: one uses a custom
rendering program and the
other uses Virtual Reality Modeling Language (VRML).  In this chapter, we will
discuses the design and implementation of our custom program.  The design and
implementation of VRML terrain rendering will be discussed in
chapter~\ref{VRML_terrain}.

At the beginning of the web technology explosion, the end user's modem speed was
slow but the general consensus was that eventually high speed networks would be
made
available to the end users.  Assuming the end user has access to a high speed
network and a parallel computer's cost effective computing power, then terrain
rendering is one of the CPU intensive applications that could be provided.  To
capitalize on these features, a terrain rendering program for a client-server
environment is needed.

We first provide an overview of a number of computer graphics techniques and the
details of our parallel rendering program.  Finally, we present our findings and
the improvements possible for our program.

\section{Computer Graphics Techniques}
\label{CG_techniques}

\begin{figure}
\begin{center}
\begin{picture}(317,274)
\psfig{figure=scan-convert.eps,width=317pt,height=274pt}
\end{picture}
\end{center}
\caption{The four sided polygon is defined by the vertices.  The order of
the vertex in the polygon's vertex list is indicated by the number.
The dash lines are the scanlines. 
\label{scan-convert}}
\end{figure}

The quality of the rendered image depends on the computer graphics
techniques used during rasterization.  Our parallel terrain rendering
program uses a scan conversion algorithm~\cite{Foley} \cite{PERogers} to
rasterize the polygons.  Our implementation of the scan conversion algorithm
assumes the polygon vertices are listed in clockwise order.
The algorithm first finds the starting scanline (scanline A).
In figure~\ref{scan-convert}, the scanlines are labeled from bottom to top,
so scanline A is at $\lceil y \rceil - 0.5$, where $y$ is the largest Y
coordinate values of the vertices and $\lceil y \rceil$ is the next integer
value larger than $y$.  The left edge (edge defined by vertices 4 and 3) and
the right edge (vertices 4 and 1) of the polygon are located on the starting
scanline.  All the pixels on the scanline between the left edge and the right
edge are visited, then the scanline is moved downward.  This process repeats
until a scanline lands on or just past a vertex, such as scanline B and C, then
a new vertex is located and its Y coordinate is used as the terminating
condition for the iteration.  In the case of scanline B, the vertex 3 is
connected to vertex 4 and 2, but only vertex 2 have a Y coordinate lower than
vertex 3; so Y coordinate of vertex 2 is used as the termination condition.
This process repeats until the left edge
intercepts the right edge (scanline D).  Then all the pixels bounded by the
polygon have been visited.  The scan conversion algorithm is a simple and
efficient method for rasterizing polygons.  In addition, it allows easy
adoption of a number of computer graphics techniques for computing the pixel
values.

\begin{figure}
\begin{center}
\begin{picture}(238,239)
\psfig{figure=gouraud-shading.eps,width=238pt,height=239pt}
\end{picture}
\end{center}
\caption{A Gouraud shaded triangle is shown.  The color of the pixel is
indicated by the circles with different shades of gray.
\label{Gouraud-shading}}
\end{figure}

\begin{figure}
\begin{center}
\begin{picture}(238,239)
\psfig{figure=phong-shading.eps,width=238pt,height=239pt}
\end{picture}
\end{center}
\caption{The arrows indicate the normal used by each pixel to compute its
color for Phong shading.  \label{Phong-shading}}
\end{figure}

Four of the computer graphics techniques, Gouraud shading, Phong shading,
bump mapping and texture mapping~\cite{Foley} \cite{PERogers},  are applicable
to terrain rendering and these techniques can also be used in conjunction with
the scan conversion algorithm.
Gouraud shading is also known as intensity interpolation shading.
The color of the object's vertices is first computed as described in step four
of the standard rendering pipeline in section~\ref{Intro_CG}.
Then the color of the pixel is obtained by linear interpolation using
the color values of the nearby vertices as shown in
figure~\ref{Gouraud-shading}.
Phong shading is normal-vector interpolation shading.
Instead of interpolating the color value of the vertices, the normal for the
pixel is calculated by interpolation using the normal of the nearby vertices
(figure~\ref{Phong-shading}).  After the
pixel's normal is obtained, the color of the pixel is computed the same way as a
vertex in Gouraud shading.  Phong shading is obviously more computational
expensive than Gouraud shading.  However, Phong shading gives more
realistically shaded rendering.
If a light source causes a bright spot to appear on an object and the bright
spot falls between the vertices, then
the bright spot will be missing from the Gouraud shaded rendering of the object.
As a result, a rotating Gouraud shaded object will have a blinking bright spot,
if it alternately falls on and between vertices.
However, for Phong shaded objects the highlight will remain unchanged.
Hence, Phong shading is well suited for specular surfaces, such as ocean
surface viewing from planetary orbit (figure~\ref{sunrise}).
Since land is not specular, Gouraud shading can be used without affecting the
image quality (figure~\ref{landmarks}).

Rendering a realistic surface texture (all the grooves, bumps and colors) is
very computationally intensive because of the large amount of information of
complex surface texture that has to be processed.  By using bump and texture
mapping, rendering effort can be reduced.  Bump mapping, which we have suggested
using in section~\ref{earth-bump-map}, can be consider as an extension of the
Phong shading, since a slight perturbation to the computed normal is added.
We recall that the object is formed by a collection of polygons, and each
polygon emulates a smooth planar surface.
The pixel normal of an Phong shaded object is obtained from the interpolation
of the polygon's vertex normals.  As a result, the normal varies uniformly
across the surface.  However, if a perturbation is added to this normal
before computing the highlight, then the color intensity will
varies accordingly.  Hence, an illusion of a bumpy surface is created.

\begin{figure}
\begin{center}
\begin{picture}(311,179)
\psfig{figure=raster-tmap.eps,width=311pt,height=179pt}
\end{picture}
\end{center}
\caption{The polygon is transformed to the screen coordinates to find out the
size and shape of it occupied section of the screen.  Then each screen pixel
within the occupied region is
backward mapped onto the polygon.  The color of the screen pixel is computed
using the color values of the texture pixel(s) under the backward mapped screen
pixel.  \label{raster-tmap}}
\end{figure}

\begin{figure}
\begin{center}
\begin{picture}(309,202)
\psfig{figure=raster_one2m.eps,width=309pt,height=202pt}
\end{picture}
\end{center}
\caption{The polygon is translated some distance closer to the viewpoint.
Each screen pixel within the occupied region is backward mapped onto the
polygon, and the area covered by the backward mapped screen pixel is smaller
than a texture pixel.\label{texture-one2m}}
\end{figure}

\begin{figure}
\begin{center}
\begin{picture}(309,183)
\psfig{figure=raster_m2one.eps,width=309pt,height=183pt}
\end{picture}
\end{center}
\caption{The polygon is rotated $45^{\circ}$ and translated some distance away
from the viewpoint. Each screen pixel within the occupied region is backward
mapped onto the polygon, and the backward mapped screen pixel area covers a
number texture pixels.\label{texture-m2one}}
\end{figure}

The use of texture mapping, a well known computer graphics technique, allows a
polygon to be decorated with an image rather than using simple interpolation of
the vertex color.  Because each pixel in the texture can have a different color,
a texture mapped polygon have higher visual complexity.
Rasterizing a texture mapped polygon with a scan conversion algorithm requires
mapping the screen pixel back to the image coordinate as show in
figure~\ref{raster-tmap}.  Depending on the distance and
orientation of the polygon, a screen pixel can be mapped onto many texture
pixels (figure~\ref{texture-one2m}) or many screen pixel can be mapped
onto one texture pixel (figure~\ref{texture-m2one}).  In additional to
the comparatively high overhead in the backward mapping, the non-unique map
also introduces the typical difficulties of the proper handling of up and down
sampling of digital signals~\cite{Wolberg-Warping}.  However, there
are a number of approximation techniques available~\cite{Watt-Animation}.
The use of texture mapped polygons trades higher
rasterization overhead for an increase in the visual complexity of the
virtual world without significant increase in the model complexity.
In the next section, we will describe a terrain rendering system that uses
some of these techniques.

\section{Terrain Rendering}

\begin{figure}
\begin{center}
\begin{picture}(377,254)
\psfig{figure=3D_arch.eps,width=377pt,height=254pt}
\end{picture}
\end{center}
\caption{The architecture of the terrain rendering system. \label{3D_arch}}
\end{figure}

The terrain rendering system is divided into three major components, the
preprocessing, the client-server program and the rendering engine, as shown in
figure~\ref{3D_arch}.  The preprocessing not only repackages the terrain data to
make it more modular, it also capitalizes on the viewing characteristics to
reduce the runtime rendering work.  The client-server program provides the user
interface for the terrain rendering system.  In this design, the server and
the rendering engine are separate processes.  Hence, the inner workings of the
rendering engine are isolated from the server.  As a result, the implementation
of the rendering engine is flexible. The parallel
computer is used not only because of the enormous size of the terrain data, but
we also want to exploit the parallelism of the terrain data and the
rendering algorithm to speed up the computation required.  In the following
sections, we will describe each of the three components in detail.

\subsection{Preprocessing}
\label{parallel-preprocessing}

The preprocessing process converts the 2D terrain data into persistent terrain
computer graphics objects, which in turn are used by the rendering
engine.  These persistent objects are generated in two steps. First the data
is repackaged into polygon-based
computer graphics objects.  Second, the color of the vertices is computed
according to a specified lighting condition.  Because the two tasks are
performed on the entire database, it is a good candidate for parallelization.

%\begin{figure}
%\begin{center}
%\begin{picture}(195,216)
%\psfig{figure=terrain-data.eps,width=195pt,height=216pt}
%\end{picture}
%\end{center}
%\caption{The format of the DEM and Landsat images \label{terrain-dataformat}}
%\end{figure}
%The structure of the data files is shown in figure~\ref{terrain-dataformat}.

\begin{figure}
\begin{center}
\begin{picture}(189,189)
\psfig{figure=100to30DEM.eps,width=189pt,height=189pt}
\end{picture}
\end{center}
\caption{Converting a 100m DEM to a 30m DEM.  The gray scale of vertices
indicates the height value.  \label{100to30DEM}}
\end{figure}

The data for our terrain rendering program is composed from two datasets, the 
Digital Elevation Model (DEM) and the Landsat sensor values.  Both datasets 
are coregistered to eliminate the visual distortion from any mismatch of the 
datasets.  The header
contains the UTM coordinates of the area and the size of data.  The
body is a 2 dimension array of short integer (2 bytes) values for the height
field and one byte for sensor image data.  The Landsat
image is divided into 7 color bands and stored in 7 datafiles.  Six of the
seven bands are used to recreate the true color image. Because we cannot find
any information on how to convert the Landsat data into the RGB color 
value of the terrain data, the sensor values pass through an ad hoc mapping
function.  This mapping function is formulated using a
trial-and-error process based on
suggestions from people at NASA Jet Propulsion Laboratory (JPL), who helped 
process the Landsat data for the JPL terrain rendering
program~\cite{Li-1992-Rendering}, and a ``looks good'' criterion.
The original DEM from United States Geological Survey (USGS) has 100m
resolution.  Because the 
Landsat data has a 30m resolution, a 30m resolution DEM was constructed
by replicating every data point in the 100m DEM two times as shown in
figure~\ref{100to30DEM}, so 
that there is a one to one mapping of the height value to the sensor values.

\begin{figure}
\begin{center}
\begin{tabular}{|c|c|c|} \hline
Input file & Size & Type \\ \hline
Satellite Image (each band) & 10.6Mbytes & Binary \\ \hline
Digital Elevation  & 21.2Mbytes & Binary \\ \hline
\end{tabular}
\end{center}
\caption{The file size of one scene satellite image and its corresponding DEM
\label{inputdata-size}}
\end{figure}

The generic rendering process calls for the transformation of the entire data
set before the active region can be found by the visible calculation.  Because
of the enormous terrain data size, the transformation will require supercomputer
resources.  The size of
the input files for one scene is given in table~\ref{inputdata-size}.  The size
of the 2D data for the entire New York state is approximately 1Gbytes.  The 
transformation of a scene of 100km $\times$ 100Km (3333 $\times$ 3333 
vertices) terrain dataset requires 133M multiplications and 133M additions.  
The entire New York state dataset has approximately 12 scenes.
Hence, real-time rendering by using the generic rendering pipeline of a 
large terrain dataset will be prohibitively expensive (at least 7.8 Gflops
typically required for 30 frames per second).

\begin{figure}
\begin{center}
\begin{picture}(414,410)
\psfig{figure=repackage.eps,width=414pt,height=410pt}
\end{picture}
\end{center}
\caption{Repackaging the terrain data set. \label{repackage}}
\end{figure}

Even though the terrain data as a whole offers fine grain parallelism only, we
can eliminate the data dependency by repackaging the entire terrain data set
into $M$ independent computer graphics objects. The 2D terrain data is divided
into a collection of 2D sub-blocks with the boundary overlapped, which
conforms to the suggestion in section~\ref{model-joining} on joining the
terrain models, as shown in figure~\ref{repackage}. 
Each sub-block is independent from each other as far as the rendering process
is concerned because there are no links crossing the section boundary.
Consequently, changing resolution or eliminating an object will not introduce
any inconsistency, such as incomplete polygons or undefined vertices.
However, these objects remain connected in the rendered image because of the
replicated data-points.  This enable an sequential rendering engine to manipulate
the virtual world in $O(M)$ instead of $O(N^{2})$, where $M$ is the total
number of terrain graphical objects and $N$ is the size of the entire terrain data.
Although the amount of work on processing of all $M$ objects is equal to 
the processing of $N^{2}$ terrain data, an algorithmic improvement can be
realized, if the images can be rendered without processing all $M$ objects.

\begin{figure}
\begin{center}
\begin{picture}(353,261)
\psfig{figure=vertices.eps,width=353pt,height=261pt}
\end{picture}
\end{center}
\caption{A $n \times m$ data is block read by one of the processing nodes.
The numbers show the index associate with that vertex. \label{vertices}}
\end{figure}

\begin{figure}
\begin{tabbing}
123\=123\=123\=123\=123\=123\=123\=123\= \kill
\>Vertex0 = 0; \\
\>index = 0; \\
\>for (i = 0; i $<$ m - 1; i++) \\
\>\>Edge[index] = (Vertex0, Vertex0 + (n - 1) + 1); \\
\>\>index++; \\
\>\>for (j = 0; j $<$ n - 1; j++) \\
\>\>\>Edge[index] = (Vertex0, Vertex0 + 1); \\
\>\>\>index++; \\
\>\>\>Edge[index] = (Vertex0, Vertex0 + (n - 1) + 2); \\
\>\>\>index++; \\
\>\>\>Edge[index] = (Vertex0 + 1, Vertex0 + (n - 1) + 2); \\
\>\>\>index++; \\
\>\>\>Vertex0++; \\
\>\>Vertex0++; \\
\>for (j = 0; j $<$ n - 1; j++) \\
\>\>Edge[index] = (Vertex0, Vertex0 + 1); \\
\>\>index++; \\
\>\>Vertex0++; \\
\end{tabbing}
\caption{The algorithm used to generate a edge list.  \label{generate-edges}}
\end{figure}

\begin{figure}
\begin{center}
\begin{picture}(329,315)
\psfig{figure=edges.eps,width=329pt,height=315pt}
\end{picture}
\end{center}
\caption{Edges are generated for $n \times m$ vertex list.
The number is the index associated with that edge. \label{edges}}
\end{figure}

\begin{figure}
\begin{tabbing}
123\=123\=123\=123\=123\=123\=123\=123\= \kill
\>Vertex0 = 0; \\
\>Edge0++; \\
\>index = 0; \\
\>for (i = 0; i $<$ m - 1; i++) \\
\>\>for (j = 0; j $<$ n - 1; j++) \\
\>\>\>Polygon[index] = \{(Vertex0, Vertex0 + (n - 1) + 1, Vertex0 + (n - 1) + 2), \\
\>\>\>\>(Edge0, Edge0 + ((n - 1) * 3) + 2, Edge0 + 2)\} \\
\>\>\>index++; \\
\>\>\>Polygon[index] = \{(Vertex0, Vertex0 + (n - 1) + 2, Vertex0 + 1), \\
\>\>\>\>(Edge0 + 2, Edge0 + 3, Edge0 + 1)\} \\
\>\>\>index++; \\
\>\>\>Vertex0++; \\
\>\>\>Edge0 += 3; \\
\>\>Vertex0++; \\
\>\>Edge0++; \\
\>for (j = 0; j $<$ n - 1; j++) \\
\>\>Polygon[index] = \{(Vertex0, Vertex0 + (n - 1) + 1, Vertex0 + (n - 1) + 2), \\
\>\>\>\>(Edge0, Edge0 + ((n - 1) * 3) + 1 + j, Edge0 + 2)\} \\
\>\>index++; \\
\>\>Polygon[index] = \{(Vertex0, Vertex0 + (n - 1) + 2, Vertex0 + 1), \\
\>\>\>\>(Edge0 + 2, Edge0 + 3, Edge0 + 1)\} \\
\>\>index++; \\
\>\>Vertex0++; \\
\>\>Edge0 += 3; \\
\end{tabbing}
\caption{The algorithm used to generate a polygon list.
\label{generate-polygons}}
\end{figure}

\begin{figure}
\begin{center}
\begin{picture}(297,297)
\psfig{figure=polygons.eps,width=297pt,height=297pt}
\end{picture}
\end{center}
\caption{Polygon are generated for the terrain computer graphics object.
The number is the index associate with that polygon. \label{polygons}}
\end{figure}

The repackaging starts by each processing node reading a $n \times m$
block, where $n$ and $m$ specify the dimension of the subsection
of the terrain data from the data files by row.  An unambiguous order is
imposed onto the points and creates the vertex list of the terrain object
(figure~\ref{vertices}).  The edges are generated using the vertices from the
neighboring column and row, and their relationship is given by the algorithm
shown in figure~\ref{generate-edges}.
Then the edges are given a unique index according to the order of creation as
shown in figure~\ref{edges}.  The polygons use the edge indices to reference
their edges. The vertex and edge list of a polygon is constructed using the
algorithm in figure~\ref{generate-polygons}.  The ordering of the polygons is
shown in figure~\ref{polygons}.
Since each processing node is reading the terrain data by the block, we 
use this opportunity to obtain the global block-block distribution
needed by the lighting calculation, which computes the color of the vertices.

\begin{figure}
\begin{center}
\begin{picture}(217,131)
\psfig{figure=lighting-angle.eps,width=217pt,height=131pt}
\end{picture}
\end{center}
\caption{The ideal lighting angle is between the two $45^{\circ}$ arc periods
\label{lighting-angle}}
\end{figure}

The viewing characteristics of terrain data are exploited, so that part of
the lighting calculation can be done in the preprocessing step and the rest
is completed at runtime.  By dividing lighting calculation into two part,
the amount of work for rendering is reduced, while the terrain rendering
program still retains the flexibility of adjusting the lighting intensity at
runtime.  The precomputed lighting is possible due to the invariance
of the terrain dataset and the lighting source (the sun) over the time of the
rendering simulation, which allows us to choose in advance the best viewing
light condition.   The terrain features are seldom 
changed within the human time scale. For example, the USGS DEM is obtained 
from digitizing existing maps that are at least 20 years old.   Some of the 
surveys for the maps were done at the beginning of the 20th century.  Because of
the satellite orbit and natural obscuring features, such as cloud cover and fog,
a clear satellite image is a chance encounter.  Hence, the best Landsat scene is
collected over the years.  This fosters a perfect opportunity to do as much work
in the preprocessing phase as possible because the preprocessing cost is spread
over the long active lifetime of the dataset.

Pre-selecting the sun position for the best viewing condition also enables
the lighting of the terrain to be computed in the preprocessing phase.
Assuming the user views the terrain on a dry and clear summer day with 12 hours
daylight, and 
the sun sweeps across the sky from east to west in a $180^{\circ}$ arc.  The sun
can only provide the best viewing light conditions in a limited range of angles.
For the 2 hour period around dawn and dusk ($30^{\circ}$ above the 
ground), the shadows of the terrain are elongated. The shading technique, which 
is used for rendering the terrain data, does not provide a good approximation 
for the elongated shadow (see section~\ref{earth-bump-map} for a more detailed
discussion).  For 
the 2 hours around noon ($\pm 15^{\circ}$ of the vertical), there is not enough 
contrast to allow the user to easily discern features in three dimensions.
For example, without any shadow, a hill to the east of the viewer looks exactly
like the hill to the west.  Hence, the user is easily lost within the virtual world
because of the lack of a quick navigational reference.  These three
time periods (total of $90^{\circ}$ arc interval) are not suitable for viewing
the terrain data.  In the remaining two $45^{\circ}$ arc periods shown in
figure~\ref{lighting-angle}, the sun is either angled slightly to the east or to
the west.  The rendered images under both lighting angle periods are similar
except that the shadow is pointing in a different direction.  Therefore, the
global light condition can be pre-selected from one of the two intervals to
provide the best viewing condition, and the most CPU intensive part of the
lighting calculation can be completed at pre-processing phase.

\begin{figure}
\begin{center}
\begin{picture}(333,333)
\psfig{figure=object-comm.eps,width=333pt,height=333pt}
\end{picture}
\end{center}
\caption{Communication pattern to share the partial normal value with the
neighboring object.  \label{object-comm}}
\end{figure}

\begin{figure}
\begin{center}
\begin{picture}(432,432)
\psfig{figure=block-object-comm.eps,width=432pt,height=432pt}
\end{picture}
\end{center}
\caption{Communication pattern to share the partial normal value with the
neighboring object with block-block data distribution.
\label{block-object-comm}}
\end{figure}

\begin{figure}
\begin{center}
\begin{picture}(460,146)
\psfig{figure=vertex-normal.eps,width=460pt,height=146pt}
\end{picture}
\end{center}
\caption{A vertex normal is the normalized sum of the adjacent polygon normals.
\label{vertex-normal}}
\end{figure}

Since terrain is non-reflective (ie, $k_{d}$ in
formula~\ref{full-light-intensity-formula} is a constant), the lighting
calculation consists of determining if the vertex is facing the light source or
not.  The light intensity value of each vertex is given by the following formula
\begin{displaymath}
\begin{array}{lcl}
I & = & I_{P} N'_{HF} \cdot L \\
\end{array}, \label{terrain-light-intensity-formula}
\end{displaymath} where $I_{P}$ is the maximum intensity value,
$N'_{HF}$ is the surface normal and $L$ is the vector pointing toward the
light source.  If the light
normal and the vertex normal oppose each other, then the vertex is lighted.
Otherwise, the vertex is dark.  The lighted or shaded condition can be
determined by the sign of the dot product of the light normal and the vertex
normal.  The direction of a vertex normal is computed by summing all the
normals of the polygons which reference the vertex
(figure~\ref{vertex-normal}).  Then the resulting vector is normalized and
assigned to the vertex.  The polygon normal is computed by using the three
vertices to form two vectors on the surface of the polygon.  The cross
product of these two vectors is the polygon normal.
Since the vertices on the boundary of an object have only a partial
sum of these direction vectors, information has to be exchanged with the
neighboring objects.  The communication pattern for exchanging information
between objects is shown in figure~\ref{object-comm}.  In order to ensure the
partial value of the vertices which exist in four objects is shared, the
communication is divided into two stages, north-south and east-west.
With block-block
distribution of the terrain data over multiple nodes, some of the neighboring
objects are located in the same node, so sharing can be done by memory copy,
thus reducing the overall communication cost.  While boundary
values require message passing between the nodes, and the
communication pattern become as shown in figure~\ref{block-object-comm}.

After the lighting condition is computed, each terrain graphics object is
written to a single file.  The object is completely encoded in three lists:
a vertex list, an edge list and a polygon list.  Each entry of the vertex
list contains the 3D coordinates and the RGB color of the vertex.  An edge list
entry has the starting and ending vertex indices.  A polygon list entry
consists of the six indices of its vertices and edges.  The structure of the
data file is shown in figure~\ref{object-file}.
After all the objects are saved, the rendering engine can use the data to
generate images.  The control of the rendering engine is provided by the
client-server program, which is describe in the next section.

\begin{figure}
\begin{center}
\begin{picture}(112,430)
\psfig{figure=obj-file.eps,width=112pt,height=430pt}
\end{picture}
\end{center}
\caption{The file structure for a persistent terrain computer graphics
object.  \label{object-file}}
\end{figure}

\subsection{Client-Server Program}
\label{client-server}

Our client-server design provides an efficient and flexible
interface for the terrain rendering engine.   The three tier arrangement allows
the encapsulation of specific tasks into independent processes, so that the
rendering process can be streamlined.   In this section, we will first describe
the client process, then the server process.

\begin{figure}
\begin{center}
\begin{picture}(432,479)
\psfig{figure=client_screen.ps,width=432pt,height=479pt}
\end{picture}
\end{center}
\caption{The user interface of bi-plane flight simulator\label{tcl-client}}
\end{figure}

\begin{figure}
\begin{center}
\begin{picture}(432,259)
\psfig{figure=view_info.eps,width=432pt,height=259pt}
\end{picture}
\end{center}
\caption{The relationship of the data entries in the {\tt view\_info} structure.
\label{view_info}}
\end{figure}

The client process provides a user interface and the application context
for the terrain rendering engine.  Possible applications include, but are not
limited to, taking a picture of the area from a specified point in space,
touring the area
with teleportation, or flying an aircraft.  However, before the rendering engine
becomes operational, the client must provide the rendering engine with the size
of the output image.  Our client is a bi-plane flight simulator (the plane type
was selected based on overall performance of the system and the smoothness of
the presentation) as shown in figure~\ref{tcl-client}.
It provides the basic motion control, such as speed,  ascent rate, descent rate,
left and right turn, for the user's virtual bi-plane and the position
of the virtual pilot's head, for example looking left, right, up and down. Our
client is implemented by using Tcl/Tk~\cite{TclTk-toolkit}
\cite{TclTk-programming} and C.  All motion controls are
achieved by adjusting the appropriate sliders.   These adjustments represent
incremental corrections to the current position and orientation of the viewer.
The information of the next view point is computed by combining the inputs
from the controls with the
previous coordinate, which is stored in {\tt view\_info} and
{\tt View\_matrix}.  The {\tt view\_info} data structure contains the following:

\begin{tabbing}
123\={\tt offset\_angle1, offset\_angle2}\=1\=123\=123\=123\=123\=123\= \kill
\>{\tt start\_x, start\_y} \> \>the anchor coordinates of the view volume 
projection\\
\>{\tt angle} \> \> the direction the viewer is facing \\
\>{\tt offset\_angle1, offset\_angle2} \> \> the width of the view volume
wedges.
\end{tabbing}

The relationship of these information values is shown in figure~\ref{view_info}.
The {\tt View\_matrix} is $4 \times 4$ $U$ matrix, which is a combined
transformation matrix described in section~\ref{Intro_CG}.
Since the rendering engine only accepts {\tt view\_info} and {\tt View\_matrix}
as inputs, it is completely unaware of how the two data structures come into
existence.  For example, the client can be converting a clickable map coordinate
into the virtual world viewing information.  In that case, the successive
coordinates that are received by the rendering engine appear random.  Or, the
client can be an aircraft simulator; hence, the successive coordinates are
incremental adjustments of the previous coordinates that are based on the user
inputs.  The terrain rendering engine can be used in both systems without
modification.

After all needed information is gathered, the data in the {\tt view\_info} and
{\tt View\_matrix} structures are sent to the server via a TCP/IP connection.
When the rendering of the image is completed, the client receives the resulting
image from the same TCP/IP connection.
To reduce the jittering of displaying successive images, our bi-plane simulator
maintains two image buffers.  While displaying one buffer, the other is
receiving a new image.  The act of receiving the image involves moving a
large block of data from the main memory to the screen buffer and doing
color quantization.  If the client has graphics rendering hardware, then the
data movement and quantization can be handled by the specialized hardware.
Otherwise, the work is done by the client software.
Furthermore, the client also needs to compute the parameters for
the future image to keep the rendering engine busy.  This requires the client to
have multi-thread capabilities. Since Tcl/Tk does not support a multi-thread
programming environment, we implemented a C widget which alternates
between reading the user input and the image from the server.
Because of the simplistic scheduling, the C widget requires proper tuning to
match the ``apparent'' rendering speed, which is affected by the combination of
network bandwidth, network traffic, the rendering hardware of the rendering
engine and the graphics display capabilities of the client.  This can be
eliminated if the client is implemented in
Java~\cite{Cornell-Core-Java} (see section~\ref{java-client} for more details).

The rendering process for the rendering engine can be divided into three stages:
\begin{enumerate}
\item receive view point from client
\item render the image
\item send image to client
\end{enumerate}
The communication and the rendering work can proceed in parallel for rendering
multiple images.  An efficient rendering process needs two threads, one to
communicate with the client and the other for rendering.  However, none of the
parallel operating systems support threads. We decided to emulate the two
threads with two UNIX processes, a server process and a rendering engine
process.

The server acts as an information broker between the client and the rendering
engine during software development and runtime.  During the development of the
client and the rendering engine software, the existence of a server provides a
logical isolation of the programming details, such
as the changes in Application Programming Interface (API), the differences of
the computer architecture and connectivity, and the limitations of the
programming
environments.  For example, Parallel Operating Environment (POE) on the IBM
SP2 computer enables the parallel program to communicate using either the high
speed network or a regular TCP/IP connection without recompilation.
A server process is used to avoid any conflicts between the
TCP/IP connection to the client and the internal communication channels between
the nodes.  The two processes communicate by using UNIX's interprocess
communication facilities.   However, POE has an
undocumented limitation, which disables most UNIX interprocess communication
facilities, except the shared memory facility.  Although this does not have a
direct effect on the client API, it has some implications on how and when the
client data is received and the image sent.  As a result, the isolation of the
client and the rendering engine reduces code complexity and runtime overhead
of the client and the rendering engine.
At runtime, the server handles the set-up and tear-down of the communication
channels, handles the sending and receiving of information, and converts the
data between network format and architecture dependent format.
With the server servicing the communication channels, the communication
overhead between the client and the rendering engine is reduced.

The three tier design of our client-server provides the rendering engine with a
programming and operating environment that is independent from the specifics
of the application.  As a result, the efficiency and portability of the
rendering engine is increased.  In the following section, we will discuss the
detail of our implementation of the rendering engine.

\subsection{Rendering}
\label{parallel-rendering}

In order to have fast rendering of a large terrain data set, a parallel
rendering engine is implemented using the MPICH-A~\cite{MPICH}, which is the
product of a join project between Argonne National Laboratory and Mississippi
State University on implementing the MPI on top of the exists machine specific
communication libraries.
The parallelism available for the terrain rendering is
in the form of object parallelism, polygon parallelism, pixel parallelism and
frame parallelism.
The object parallelism derives from the characteristics of the repackaged
terrain object, such
as non-overlapping, having an unambiguous depth order and modularity.
The characteristics of non-overlapping and having an unambiguous depth order
enable the arbitrary rendering order of the objects, while ensuring the same
image is generated.  Their modularity allows adjacent terrain objects be
distributed and rendered without communication. 
Each terrain computer graphics object consists of hundreds of polygons.
Although these polygons share boundary vertices, we can uncouple the dependency
by replicating the shared vertices.  This allows us to treat the polygons as a
collection of smaller independent computer graphics objects.
As in all other graphics applications, all pixels in an image are totally
independent from each other.  Hence, the pixels can be rendered in arbitrary
order or in parallel.   Another source of parallelism exists in time, when the
terrain engine is used to generate multiple images with correlated viewpoints
on the same data.  For example, a flight simulation program has to generate
multiple images per second and the viewpoints of the images are aligned along
the flight path.  If the aircraft does not constantly changing its flight
path, then the terrain rendering program can predict the
next few viewpoints and render the images for different frames concurrently.
Our current implementation of parallel terrain rendering engine can exploit
object and pixel parallelism, but not the time-slice parallelism.

\begin{figure}
\begin{center}
\begin{picture}(230,432)
\psfig{figure=rendering_engine.eps,width=230pt,height=432pt}
\end{picture}
\end{center}
\caption{The pipeline of the terrain rendering.
\label{rendering-engine}} \end{figure}

The parallel terrain rendering engine generates each image by using the user
specified viewpoint and a distributed terrain object database.  First, the
viewpoint is boardcast to all processors by using {\tt MPI\_bcast}.
Then the needed objects are selected from all the terrain objects. These objects
are concurrently rendered onto a distributed image buffer by a collection of
sequential rendering pipelines (figure~\ref{rendering-engine}) to exploit the
inherent object parallelism.  Since the lighting calculation is done in
the preprocessing, the sequential rendering pipeline is a shortened version of
the standard rendering pipeline described in section~\ref{Intro_CG}.  Using
the standard rendering algorithm allows
the pipeline to exploit computer graphics acceleration hardware, if it exists.

\begin{table}
\begin{center}
\begin{tabular}{|c|c|c|} \hline
Runtime file & Size & Type \\ \hline
One scene & 2.32Gbytes & Binary \\ \hline
New York state & 27.9Gbytes & Binary \\ \hline
\end{tabular}
\end{center}
\caption{The estimated size of the persistent terrain computer graphics
objects files. \label{persistent-size}}
\end{table}

As in the first step of a standard pipeline, our parallel rendering pipeline
first selects a list of objects from the distributed 2D array of computer
graphics objects (the virtual world).  Distributing the 2D array is essential
because the file size of the persistent terrain computer graphics objects
increases by a factor of approximately 15 compared to the original
2D terrain data.  The estimated sizes of the runtime files are given
in table~\ref{persistent-size}. The increase in data size comes from the
explicit encoding of the coordinate
of the height values in the 3D world, and the encoding of the edges and
polygons.  Although some of these edge and polygon information can be generated
at runtime, they are explicitly written to the files, so that we
can estimate the aggregate memory used by the terrain computer graphics objects.
Although not all the data is needed at the same time, the data is still loaded
in core to avoid slow disk access or the extra work of regenerating the
information as it is done by the preprocessing program.  As a result, only a
parallel computer offers enough aggregate memory and CPU power to render the
image using the full dataset. Even then, a coarse visible
calculation is essential to reduce the number of objects that need to be
process.

The coarse visible calculation involves
rasterizing a 2D projection of the view volume (bounding polygon of all
visible data) onto the terrain data grid.  All the terrain objects that are
fully or partially within the bounding polygon are assumed to be visible.
If the viewer orientation is limited as specified in
section~\ref{visible-calculation},
the 2D view volume projection onto the data plane is always triangular.
Because of the limited size of the terrain data grid, the bounding edges of the
grid may change the original triangular polygons, so that they acquire a larger
number of edges.  There are a total of three possible polygon shapes without
considering the four-fold rotation
symmetry as shown in figure~\ref{coarse-visible-cases}.  The limited number of
configurations allows us to hardwire the rasterization to speed up the
computation.   The algorithm of the coarse visibility calculation is shown in
figure~\ref{coarse-visible-alg}.  After the shape of the view volume projection
is determined, a modified scan conversion algorithm
(section~\ref{CG_techniques}) is used to mark all the ``visible'' terrain
objects.

\begin{figure}
\begin{center}
\begin{picture}(432,125)
\psfig{figure=coarse-vis-cases.eps,width=432pt,height=125pt}
\end{picture}
\end{center}
\caption{The three possible cases of coarse visibility calculation is shown.
The gray area is the visible region.
\label{coarse-visible-cases}}
\end{figure}

\begin{figure}
\begin{tabbing}
123\=123\=123\=123\=123\=123\=123\=123\= \kill
\>Compute Line 1 intersection point with the boundary. \\
\>Compute Line 2 intersection point with the boundary. \\
\>Mark grid sections near intersection point of Line 1 and 2 visible. \\
\>if (Line 1 and Line 2 intersect the same boundary line) \\
\>\>Visible bounding polygon is a triangle. \\
\>else if (Line 1 and Line 2 intersect adjacent boundary line) \\
\>\>Visible bounding polygon is a four sided polygon. \\
\>else if (Line 1 and Line 2 intersect opposite boundary line) \\
\>\>Visible bounding polygon is a five sided polygon. \\
\end{tabbing}
\caption{The coarse visible calculation algorithm.
%Line 1 and 2 is shown in figure~\ref{coarse-visible-example}
\label{coarse-visible-alg}}
\end{figure}

\begin{figure}
\begin{center}
\begin{picture}(427,393)
\psfig{figure=coarse-vis-ex.eps,width=427pt,height=393pt}
\end{picture}
\end{center}
\caption{An example of the coarse visible calculation
\label{coarse-visible-example}}
\end{figure}

The coarse visible algorithm first performs three intersection calculations
using the vertices of the triangular view volume projection to determine the
foot print of the projection.  Then, it does two intersection calculations for
the left and
right boundary of the view volume on each row of the terrain objects within the
visible foot print.  In the worst case, the foot print can cross all rows of
terrain objects and cover the entire terrain data grid.  As a result, the coarse
visible algorithm does at most $2N + 3$ intersection calculations
and marks at most $N^{2}$ visible objects, where $N$ is the size of the
active terrain data grid, as shown in figure~\ref{coarse-visible-example}.
Each intersection calculation involves 11
floating point additions, 10 floating point multiplications and 1 floating
point division. The example shown in figure~\ref{coarse-visible-example}
requires 37 intersection calculations to eliminate 133 grid sections.  In this
case, the average cost of the visibility calculation per eliminated terrain
object is 6 floating point additions and 6 floating point multiplications.
Since the rendering cost of a terrain object is much more than the work done to
eliminate it, the coarse visibility calculation
is a  highly cost effective way to reduce the rendering resources required.

\begin{figure}
\begin{center}
\begin{picture}(263,219)
\psfig{figure=dataprogram.eps,width=263pt,height=219pt}
\end{picture}
\end{center}
\caption{The polygon, edge and vertex arrays of the rendering program.
\label{data-program}}
\end{figure}

\begin{figure}
\begin{center}
\begin{picture}(340,268)
\psfig{figure=datastructure.eps,width=340pt,height=268pt}
\end{picture}
\end{center}
\caption{Each polygon is composed of three edges and three vertices.
The arrow represents the link between the data structures.
\label{polygon-internal}}
\end{figure}

The selected objects are rendered by the sequential
rendering pipeline on the computational node, where the objects reside.
Computational nodes on high performance parallel computer are designed
to run numerically intensive programs, and the architecture is optimized
for matrix-vector operations.  Hence, efficient rendering requires the
terrain data to be organized in vector form, so the terrain object's polygons,
edges and vertices are stored in three 1D arrays
(figure~\ref{data-program}).  Organizing the vertices into a long vector works
well for the transformation, which is handled as multiplication of a
$4 \times 4$ matrix and $4 \times l$ matrix, where $l = n \times m$ is
the number of vertices in a terrain object, $n$ and $m$ is the dimension of the
DEM block.  However, the vectorized edge and polygon list are not as effective
because edges and polygons are higher level objects, and they reference their
vertices by links as shown in figure~\ref{polygon-internal}.  Each edge has
access to the two endpoints.  A polygon has access to its edges and vertices.
The current vertex storage and labeling algorithm in
section~\ref{parallel-preprocessing} does not consider the access pattern at the
edge and polygon level.  As shown in figure~\ref{data-program},  the vertices of
a polygon can be located far apart.
If the vertices are separated widely enough, then it may be possible to induce
a cache miss for the rendering of each polygon.  On the other hand, smaller
terrain objects incur additional overhead on the coarse visible calculation.
As a result, the optimal size of the terrain object must be turned for a
particular parallel computer.

After rendering all the local computer graphics objects, each node has a
partial image.  Each of the intermediate images consist of 2D arrays of pixel,
which have the following four components:

\begin{tabbing}
123\={\tt Red, Green, Blue}\=1\=123\=123\=123\=123\=123\= \kill
\>{\tt Red, Green, Blue} \> \>the three color components of the pixel\\
\>{\tt Z\_buffer} \> \> the depth information of the pixel\\
\end{tabbing}

These images are combined using {\tt MPI\_Reduce} using a custom
combine function that operates on the pixel.  This function first compares the
depth of the pixel.  If the
``new'' pixel is closer than the ``old'' pixel, then the color of the new pixel
is copied to the final image buffer.  The pixel parallelism is exploited, when
the multiple intermediate images are combined at the same time during the
operation of {\tt MPI\_Reduce}~\cite{MPICH}.   In the next section, we will
show details of our findings for the parallel rendering engine.

\section{Findings}
\label{finding}

In this section, we present the performance of the client-server program and the
parallel graphics rendering engine.  All of these tests
were done using an IBM SP2~\cite{IBM-SP2} with 8 narrow nodes, which communicate
through the high performance switch.  A shared Ethernet local area network
(LAN) with the maximum throughput of 10Mbits per second connects the server, one
of the SP2 nodes, with the client, a Sun IPC
workstation.  This network connection allows a theoretical peak throughput
of 1.5 frames per second for transferring uncompressed 600 by 480 24-bit RGB
(0.82 Mbytes) images.  In our tests, the ``apparent'' frame rate was
much lower, partly because the network was shared by other computers.
For comparison, a Moving Pictures Experts Group (MPEG) I
video stream requires only 1.5Mbits per second to maintain National Television
Standards Committee (NTSC) quality video
(approximately 352 by 240 resolution).  A 720 $\times$ 480 resolution 
MPEG II requires bandwidth from 4Mbits to 9Mbits per second.
As a result, our image delivery problem can be solved by
using a higher bandwidth network and by compressing the images before shipping.
Tradeoffs have to be made between the benefit of reducing the required bandwidth
and the time for the client to uncompress the image.
In section~\ref{client-server}, we speculated that the client's lack of
computer graphics hardware support might slow down the presentation of the
images.  Our tests show that the client can decode MPEG I video at the
rate of 10 to 12 frames per second.  It can also copy up to four
600 $\times$ 480
RGB images per second onto the screen buffers.  As a result, the
theoretical peak performance of the client-server program is limited by the
network bandwidth.

\begin{table}
\begin{center}
\begin{tabular}{|l|c|c|c|c|c|c|c|c|} \hline
Node & Number of objects per picture & Rendering Time (sec) & Communication Time \\ \hline
(0, 0) & 34 $\pm$ 46 & 0.28 $\pm$ 0.35 & 0.93 $\pm$ 1.9 \\ \hline
(0, 1) & 35 $\pm$ 48 & 0.27 $\pm$ 0.32 & 2.3 $\pm$ 3.43 \\ \hline
(1, 0) & 35 $\pm$ 48 & 0.27 $\pm$ 0.30 & 2.5 $\pm$ 7.2 \\ \hline
(1, 1) & 37 $\pm$ 49 & 0.29 $\pm$ 0.31 & 0.24 $\pm$ 0.44 \\ \hline
(2, 0) & 37 $\pm$ 50 & 0.29 $\pm$ 0.31 & 2.4 $\pm$ 3.9 \\ \hline
(2, 1) & 39 $\pm$ 50 & 0.30 $\pm$ 1.3 & 0.33 $\pm$ 0.73 \\ \hline
(3, 0) & 38 $\pm$ 50 & 0.31 $\pm$ 3.5 & 0.61 $\pm$ 0.65 \\ \hline
(3, 1) & 40 $\pm$ 51 & 0.33 $\pm$ 3.5 & 0.19 $\pm$ 0.03 \\ \hline
\end{tabular}
\end{center}
\caption{Performance of the individual parallel rendering pipeline.
The statistic is collected using 800 $\times$ 800 height values divided into a
32 $\times$ 32 grid.
\label{load_balance}}
\end{table}

The performance of a graphics rendering engine is typically measured by its
\begin{eqnarray}
\mbox{\small Frame Rate} & = & \frac{\mbox{\small 1}}{\mbox{\small Total time per image}}
\end{eqnarray}
and
\begin{eqnarray} 
\mbox{\small Polygon Rendering Rate} & = &
\frac{\mbox{\small Number of objects per image} \times \mbox{\small Polygons per object}}
{\mbox{\small Total time per image}}.
\end{eqnarray}
Although these measurements are often quoted for comparing commercial products,
their validity is often questioned because of the extensive impact of one or
more arbitrary testing parameters~\cite{benchmark1} \cite{benchmark2},
such as lighting model, viewpoints
and data set.  As shown in section~\ref{CG_techniques}, the computation
of different lighting models involves different amounts of calculations. In our
rendering engine, Gouraud shading is used, but some of the calculations are
done off-line, so our rendering speed-up is partly due to the reduced
workload.  The effect of the viewpoint of the rendered images is also difficult 
to qualify, yet has far reaching
impact on the benchmarking process. If the viewpoint is far from the objects,
then the polygons are represented by only a small number of screen pixels.
Since the algorithmic complexity of polygon rasterization is bounded by the
number of screen pixels covered, a small polygon is usually rendered
faster.  Furthermore, the combination of viewpoint and the local topological
features of the terrain data can affect the number of visible polygons.  Suppose
the terrain data is a 2D uniform samples of height values.  If a mountain
is centered in the data set and the viewpoint is located such that the mountain 
blocks the view to the rest of the polygons, then the rendering speed of that
image is atypically faster than rendering a plain represented by the same
number of polygons, since a reduced amount of work is required to generate the
images.  For our benchmarking of the parallel rendering engine, we used a flight
path consisting of a total 311 viewpoints, cutting diagonally
across the terrain data grid, so that the timing is done under extreme
conditions, since in some images almost all the terrain objects are visible
and in some images only a small number of them are visible. This, however, led
to the large uncertainty of the timing data.

In all our performance testing of the parallel rendering engine, the same
$800 \times 800$ points, a region of 24Km $\times$ 24Km relatively flat terrain
near Albany, New York, was used to generate the performance statistics.  The
terrain data is divided into 1024 terrain objects.  Each terrain object has
approximately 1250 polygons.  An average of 295 terrain objects, or
368,750 polygons, are rendered for each image.  The average total
time for rendering an image is the sum of the maximum average rendering time,
0.33 second, and the communication time at node $(0, 0)$, 0.93 second.  The
communication time of node $(0, 0)$ is used because the final image is
gathered at node $(0, 0)$ before serving to the client.  As a result, the
communication timing of node $(0, 0)$ shows the true communication time
needed to collect all the partial images.  The rendering speed of
node $(0, 0)$ is slightly slower because it also handles the extra work of
sending the image to the client.  The parallel rendering engine on 8
processors has achieved a stable frame rate of 0.77 frames per second for an
$800 \times 800$ points and it is capable of handling more than 280,000 polygons
per second.  By comparison, the SGI ONYX~\cite{SGI-ONYX} with a
RealityEngine II have peak rendering rate of 1.8M triangle per second.  Although
our parallel rendering engine only achieves 16\% of the peak performance of the
RealityEngine II, the software still have room for improvement (see
section~\ref{improvement}).

\begin{figure}
\begin{center}
\begin{picture}(432,127)
\psfig{figure=cyclic-cyclic.eps,width=432pt,height=127pt}
\end{picture}
\end{center}
\caption{Cyclic cyclic distribution of the terrain object \label{cyclic-cyclic}}
\end{figure}

\begin{table}
\begin{center}
\begin{tabular}{|l|c|c|c|c|} \hline
 & & \multicolumn{3}{c|}{test case} \\ \hline
Processor & Terrain object grid & 1 & 2 & 3 \\ \hline
 1 & 25 (31\%) & 6 (29\%) & 8 (29\%) & 14 (29\%) \\ \hline
 2 & 20 (25\%) & 6 (29\%) & 7 (25\%) & 12 (25\%) \\ \hline
 3 & 20 (25\%) & 4 (19\%) & 7 (25\%) & 12 (25\%) \\ \hline
 4 & 16 (20\%) & 5 (24\%) & 6 (21\%) & 10 (21\%) \\ \hline
\end{tabular}
\end{center}
\caption{Load distribution on 4 processors with cyclic-cyclic
distribution pattern. The three test cases are given in
figure~\ref{cyclic-cyclic}.  \label{cyclic-cyclic-cal}}
\end{table}

The performance of the parallel rendering engine is disappointing because of
the load imbalance of the rendering pipelines, the low performance of the
software rendering pipeline and the time spent on gathering of the images. The
load balance of the rendering pipelines depends on the distribution of the
terrain objects.  The two dimensional grid of terrain objects are statically
cyclic-cyclic distributed among the computational nodes as shown in
figure~\ref{cyclic-cyclic}.  However, during runtime only a small sections
within the view volume is being rendered, the load of the computation nodes
become non-uniform.   Although the
cyclic-cyclic distribution does not guarantee optimal distribution at runtime,
it is the best static distribution for the terrain data because the initial
viewing orientation is randomly chosen by the user.  Hence, the initial
projection of the view volume can be consider as a random selection of a set of
adjacent terrain objects.  Since cyclic-cyclic distribution minimized any
clustering, this randomly selected set of
terrain objects are less likely to be located on the same computing node.
As shown in table~\ref{cyclic-cyclic-cal}, the cyclic-cyclic distribution
pattern for the terrain
object grid in figure~\ref{cyclic-cyclic} is not optimal.  Processor 1 has 9
terrain objects more than processor 4.  However, the load imbalance of the
three test cases are almost the same(table~\ref{cyclic-cyclic-cal}).  The
imbalance is primary due to the unbalance distribution of the terrain object
grid.  The cyclic-cyclic distribution offers the best static distribution for
the distributed terrain computer object database.
In the timing tests, the terrain data was divided into a 32 $\times$ 32 terrain
grid. The 8 nodes are mapped into a 4 $\times$ 2 processor grid, so every
processor handles 128 terrain objects.  The timing data in
table~\ref{load_balance} shows that the rendering engine does not have a
perfectly balanced load across the generation of the 311 test images. The
subsequent views and the geometry of their viewing volume projections are
very similar (section~\ref{visible-calculation}) at least over short time
period, but still many frames is rendered.  This implies that the loads
of the rendering pipelines do not vary rapidly from one image to the next image.
Hence, if an imbalanced load exists, its impact is extended to multiple
images.  The communication times of the other nodes indicate that some node
consistently finishes the rendering phase early and they are
waiting to complete the communication (table~\ref{load_balance}).  Hence, load
imbalance does happen. An adaptive load balancing algorithm is needed to
speed-up the parallel rendering pipeline.

\begin{table}
\begin{center}
\begin{tabular}{|l|c|c|c|c|} \hline
Grid size & $8 \times 8$ & $12 \times 12$ & $14 \times 14$ & $16 \times 16$ \\ \hline
Rendering Time (sec) &
0.55 $\pm$ 0.23 & 0.54 $\pm$ 0.23 & 0.54 $\pm$ 0.23 & 0.54 $\pm$ 0.23 \\ \hline
Write Time (sec) &
0.043 $\pm$ 0.012 & 0.045 $\pm$ 0.012 & 0.044 $\pm$ 0.012 & 0.047 $\pm$ 0.013 \\ \hline
Total Time (sec) & 0.59 & 0.59 & 0.58 & 0.59 \\ \hline
Frame rate & 1.7 & 1.7 & 1.7 & 1.7 \\ \hline
\end{tabular}
\end{center}
\caption{Performance of the sequential rendering engine measure by frame
rate for 400 $\times$ 400 height values\label{seq_time_FR}}
\end{table}

\begin{table}
\begin{center}
\begin{tabular}{|l|c|c|c|c|} \hline
Grid size & $8 \times 8$ & $12 \times 12$ & $14 \times 14$ & $16 \times 16$ \\ \hline
Polygons per object & 5000 & 2222 & 1633 & 1250 \\ \hline
Number of objects per picture &
29 $\pm$ 13 & 60 $\pm$ 29 & 80 $\pm$ 38 & 100 $\pm$ 49 \\ \hline
Polygons per picture & 150, 000 & 130, 000 & 130, 000 & 130, 000 \\ \hline
Total Time (sec) per picture & 0.59 & 0.59 & 0.58 & 0.59 \\ \hline
Polygons per second & 250, 000 & 230, 000 & 230, 000 & 210, 000 \\ \hline 
\end{tabular}
\end{center}
\caption{Performance of the sequential rendering engine measure by polygon
rendering rate for 400 $\times$ 400 height values\label{seq_time_PR}}
\end{table}

\begin{figure}
\begin{center}
\begin{picture}(432,272)
\psfig{figure=time1.ps,width=432pt,height=272pt}
\end{picture}
\end{center}
\caption{Breakdown of time usage in sequential rendering engine
\label{seq_time_breakdown}}
\end{figure}

The performance of the sequential rendering pipeline also plays a significant
role in the performance of the parallel rendering engine.  The performance
statistics of the sequential pipeline are collected by using a flight path,
consisting of a total 114 viewpoints, cutting diagonally across the terrain data
grid of $400 \times 400$ points, a region of 12Km $\times$ 12Km relatively flat
terrain near Albany, New York.  The timing data shows that the sequential
rendering pipeline has a stable frame rate of 1.7 frame per second for
640 $\times$ 400 images across four different terrain object sizes as shown in
table~\ref{seq_time_FR}.  The pipeline is capable of handling more than
200,000 polygons per second (table~\ref{seq_time_PR}). Contrary to our
speculation in section~\ref{parallel-rendering}, the widely scatter vertices
did not induce cache
misses.  A more detailed breakdown of the rendering process is shown 
in figure~\ref{seq_time_breakdown}.  The data is normalized to eliminate the
bias cause by the {\tt gprof} performance logging functions, so that a clear
picture can be uncovered from the data.  The rendering pipeline spent
most of its effort in scan conversion of triangles, which consists of
subroutines {\tt triangle\_scan\_convert}, {\tt scanline},
{\tt incrementalize\_y}, {\tt incrementalize\_x} and {\tt increment} and
coordinate transformation, {\tt process\_vertices}.  Consequently, hardware
rendering or an assembly coded scan conversion function will dramatically
improve the performance of the rendering pipeline.

The small speed-up of the parallel rendering engine is mainly due to the
inefficient communication to assemble the final image before shipping to the
client.  The gathering of the final image is currently done by using a single
{\tt MPI\_Reduce} call on the 600 $\times$ 480
Z-buffer and RGB values (3.5Mbytes per processor).  A user defined combine
function is used to replace the pixels further away from the viewpoint
in the received image with closer-by pixels from the local image buffer.
The operation has $O(MN)$ complexity, where $M$ and $N$ is the width and
height of the image.  This communication takes 0.93 second  to
complete, which is about 75\% of the total rendering time of a image.  Better
communication performance can dramatically improve the performance of the
parallel rendering engine.  In the next section, we will show how to improve
the client-server program and the parallel rendering engine.

%(the time spent is split evenly between sending the message and processing
%the Z-buffers)

\section{Improvement}
\label{improvement}

\begin{figure}
\begin{center}
\begin{picture}(432,152)
\psfig{figure=java_client_server.eps,width=432pt,height=152pt}
\end{picture}
\end{center}
\caption{The architecture of the Java client and server.
\label{java_client_server}}
\end{figure}

The custom terrain rendering system can be dramatically improved by
reengineering all of its components.  The client-server program can readily
use the features that offer by the Java programming language.  The rendering
pipeline can benefit from improvements on load balance, communication,
rendering technique and concurrency exploitation.

\subsection{Java Client-Server Program}
\label{java-client}

As shown in section~\ref{client-server}, the implementation of the client-server
program requires multi-threading, communication, graphical user interface and 2D
image manipulation capabilities.  The Java programming language
from Sun Microsystems~\cite{Cornell-Core-Java} provides all these; in addition,
the compiled client program can be executed on different computer architectures.
The multi-threaded program architecture of the client and the server is very
similar, with a reverse flow of information.
The client uses one thread to intercept the user's input, to compute
the viewpoint and to send it to the server.  A matching thread on the server
side receives the viewpoints and sends them to the rendering engine.
The resulting images from the rendering engine are first passed to the image
handling threads on the server then to the client's counterpart.  Hence, the
client-server program operates like a conveyor of information between the user
and the rendering engine as shown in figure~\ref{java_client_server}.

Since the client and the server are loosely coupled, the information flow
between them must be self-regulated.  By using buffered blocking I/O in all
three components, the regulation of the information flow can be accomplished
automatically.  If the rendering engine's viewpoint pipeline is full and blocks
any further entry, then the server thread is suspended while writing to the
pipeline.  In turn, the client thread suspends its sampling of the inputs, when
its I/O buffer is full.  Any congestion will be automatically resolved because
only the viewpoint handling threads are suspended and the image handling
threads remain in operation.  When the rendering engine completes its work
and passes it along to the server,  the rendering engine is available to accept
a new viewpoint and  the congestion is relieved.  In addition, the use of
multi-threading and multi-image buffers enable image compression techniques
to be easily incorporated, so that the network bandwidth requirement can be
reduced.  The increased delay, caused by the compressing and decompressing
the images, is disguised as the inertia of the aircraft.

\subsection{Load Balancing Improvement}

To address the parallel rendering engine's shortcomings of load imbalance, slow
rendering and long communication time, improvements of resource management
and rendering algorithms are needed. The rendering pipelines
are considered load balanced, if they complete their work for a particular image
at the same time.   Since the rendering time of a computer graphics object can
vary according to the viewing parameters, an adaptive load balancing
algorithm, which redistributes workload according to rendering time, is needed.
A master-slave organization of the computation nodes allows the realization of
workload redistribution.  Initially, the master assigns work to the
slaves.  When a slave completes its work and returns the result, new work is
assigned by the master.  If a slave by chance consistently receives
short jobs, then it gets more jobs.  Consequently, the workload of the
pipelines is adjusted based upon the time that is required to complete it.
An additional advantage of this algorithm is that each pipeline only needs to
keep one or two terrain objects in memory through out its life time.

\begin{figure}
\begin{center}
\begin{picture}(432,204)
\psfig{figure=load_balance1.eps,width=432pt,height=204pt}
\end{picture}
\end{center}
\caption{Adaptive load balance with master-slave organization.  Terrain
objects are redistributed before generating each picture. Viewpoints differ
by turning $4^{\circ}$ to the right.  \label{load_balance1}}
\end{figure}

\begin{table}
\begin{center}
\begin{tabular}{|l|c|c|c|c|} \hline
Frame Rate &  Redistribution cost per image \\ \hline
30 & 250 \% \\ \hline
20 & 170 \% \\ \hline
10 & 83 \% \\ \hline
5 & 42 \% \\ \hline
\end{tabular}
\end{center}
\caption{Redistribution cost per image for different frame rate.
\label{rdist_cost_relative}}
\end{table}

\begin{figure}
\begin{center}
\begin{picture}(432,204)
\psfig{figure=load_balance2.eps,width=432pt,height=204pt}
\end{picture}
\end{center}
\caption{Adaptive incremental load balance of the terrain objects operates with
only minimal redistribution. Viewpoints are differ by turning
$4^{\circ}$ to the right.  \label{load_balance2}}
\end{figure}

To implement this runtime load distribution strategy, the terrain objects must
be redistributed to different pipelines as the need arises.
Figure~\ref{load_balance1} shows an instance of the
distribution of the terrain objects for two consecutive viewpoints.
The redistribution cost per image is given by the following formula,
\begin{eqnarray}
\lefteqn{\mbox{Redistribution cost per image}} \nonumber \\ 
& = & \mbox{Redistribution time per object} \times
\left\lceil\frac{\mbox{Number of objects per image}}{\mbox{number of pipelines}}
\right\rceil.
\end{eqnarray}
The 1024 terrain objects dataset from section~\ref{finding} is used to
illustrate the effectiveness of the load redistribution algorithm.  An estimated
per object redistribution time is computed by using the size of the
terrain object (73Kbytes) and the communication speed for the SP2 high
performance switch (32Mbytes per second). Hence, the estimated
adaptive load balance cost is 0.083 second per image. Since this load balance
cost is constant across different
frame rates, the relative time spent on load balancing per image in
table~\ref{rdist_cost_relative} is a better indicator for the effectiveness of
the algorithm.  This adaptive load balance algorithm is a failure because
of the extremely high overhead for redistributing the objects and because
real-time rendering is very restrictive everything must be done very quickly.
However, the
redistribution cost can be significantly reduced by exploiting the locality of
the viewing parameters.  As shown in section~\ref{finding}, the changes of
the view volume projection vary slightly from one viewpoint to the next.
Instead of resending the same terrain objects to the pipelines, the master
can inform the slaves to reuse the existing data and only distribute the new
terrain object to the slave with the lease work.  Using this
incremental  redistribution on the above dataset, each slave has to maintain
37 terrain objects in memory.  The current distributed terrain object
database, on the other hand, stores all 1024 terrain objects in core on the 8
computing nodes.  Hence, the incremental  redistribution requires substantially
less than the memory.  The slave
also records the time required to render each terrain object, and this
information is used by the master to predict the work load of the pipelines
for the next picture.  While the slaves are
busy, the master can prepare the ``new'' terrain objects and distribute the
terrain objects to the slaves.   This enables the
overlapping of the redistribution (communication) with rendering except
for the first image. An example of this adaptive incremental load balancing
strategy is shown in figure~\ref{load_balance2} for two consecutive viewpoints.
After retiring two and adding three new terrain objects, the load
of the pipelines remains balance.  Consequently, the adaptive incremental load
balance algorithm can be used to maintain a balance load for the pipelines with
minimal redistribution cost.

\subsection{Communication Improvement}

In the current rendering pipeline implementation, the communication for the
final gathering of the sub-images dominates the overall image generation
time.  A way to
reduce communication time is to send only the modified portion of the
Z-buffer along with the pixel values.  This communication pattern is natural
for the master-slaves resource organization because the slaves have to report
to the master at the completion of their work.  Since one message per terrain
object is send to the master, the overhead of sending the message , that is,
the start-up cost of sending the message and size of the message, must be
minimized.
The communication start-up cost for the master-slave organization can be
reduced by using the MPI persistent communication request.  After each
master-slave pair is connected by a persistent communication channel,
the cost of sending an additional message is only proportional to the size
of the message divided by the communication bandwidth, ie no additional
latency is incurred by the connecting the communication channel.
To further reduce the communication cost, the size of the message must also
decreased.  Instead of sending a message, which is proportional to the size of
the entire image, from every computational nodes of the rendering pipeline,
each slave sends a number of smaller messages.  For each message, the slave
first compute the size of the screen area that occupied by the terrain object,
which will
be rendered.  Then a rendering/message buffer of that size is allocated to house
the rendering result.  After the completion of the rendering, the content of
this buffer is sent to the master.  Consequently, the communication cost of
assembling the final is minimized by reducing communication overhead and 
sending minimal amount of information.

\subsection{Rendering Technique Improvement}

In section~\ref{finding}, the use of computer rendering hardware and/or
assembly code was suggested to reduce the rendering time.  However, this is
not a scalable solution because the terrain data grows as $O(\frac{1}{N^{2}})$,
where $N$ is the resolution of the satellite image.  Our current DEM and
Landsat image is 100m and 30m resolution, respectively.  During
Supercomputing 95, Peggy Lee demonstrated the JPL terrain rendering program
with a fly-by tour of San Diego using a Clay supercomputer.  The tour was
memorable because a 3m resolution Russian declassified satellite image of the
urban San Diego area was used.  If 3m resolution terrain data is used
for the entire new York state, then it is possible to overwhelm the largest
parallel machine with moderate improvement of the resolution.  An algorithmic
improvement is
needed and can be achieved by using the wavelet terrain model optimization from
chapter~\ref{multi-terrain}.  The wavelet terrain model optimization takes
advantage of the fact that the physical characteristics of the terrain data do
not depends on the resolution of our remote sensing capability.  For example,
a lake, which currently represented by thousands of sampling points
(100m$^{2}$ per pixel) is still the same lake, if it is represented by hundreds
of millions (1m$^{2}$ per pixel) of sample
points.  The difference of the two representations is in the detail of the
shoreline.  However, the current high resolution monitor can only display
a few million pixels, so only a small portion of the original data can be
displayed.  An intelligent algorithm, on the other hand, can discard
redundant information and unwanted detail.  As a result,
the terrain model of the visible and the neighboring regions (the active region)
can be modeled with less polygons while retaining the area's distinctive terrain
features. Because the model is optimized while maintaining sufficient
information to demonstrate the unique characteristic terrain features, the
optimized model can be too complex or too simplistic for the capability of the
rendering engine and the display equipment. By using texture mapped polygons
for the terrain model, this gap between model complexity and equipment
capability can be resolved by adjusting the resolution of the texture during the
rasterization (see section~\ref{CG_techniques} for the description of
rendering texture mapped polygon).  Although the rendering cost of  a
texture mapped polygon is more expensive than a Gouraud shaded
polygon, a terrain model using texture mapped polygons before wavlet
optimization have 9 times less polygons than a terrain model using Gouraud
shaded polygons. Consequently, by using wavelet terrain model optimization
and texture mapped polygons, the per image rendering time can be reduced
substantially without compromising the quality of the resulting images.  Using
the master-slave resource organization with wavelet optimization and texture
mapped polygon terrain model, the frame rate of the custom terrain rendering
system can be quite possibly doubled.

\begin{figure}
\begin{center}
\begin{picture}(432,248)
\psfig{figure=predicting.eps,width=432pt,height=248pt}
\end{picture}
\end{center}
\caption{Multiple viewpoint is rendered in parallel.  \label{predicting}}
\end{figure}

\subsection{Concurrency Improvement}

Even with all the above proposed improvements of the rendering engine are
implemented successfully, the parallel rendering engine will most likely
produce less than 30 frames per second.  The frame rate of the rendering
engine can be improved by
rendering multiple viewpoints in parallel.  The use of the two rendering
pipelines is shown in figure~\ref{predicting}.  The second viewpoint are
obtained by predicting the location of the user's craft based on its movement
history.  The two viewpoints are rendered by two identical pipelines at once,
both images should be complete more or less the same time
because of the similarity in the viewing parameters.  Then the
two images are sent to the client in the correct order.  The second image is
temporary held in a buffer, while the server is in the process of delivering the
first image.  During the rendering time, the client samples the user inputs and
uses their values to compute the next viewpoint.  This
practice can be extended to use more than two pipelines simultaneously.
Since the user inputs are only sampled at a slower rate of the image
generation, the aircraft's control will appear sluggish.  As a result, the
additional pipelines are beneficial, when the first pipeline can no longer
profit from additional rendering resources.

