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

We have presented an algorithm to make the terrain data more
manageable in chapter~\ref{multi-terrain}.  There are two ways to use these
terrain data to generate the
images for terrain visualization, one uses a custom 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 discuses in
chapter~\ref{VRML_terrain}.

At the beginning of the web technology explosion, the end user's modem speed was
slow and the general consensus is that high speed networks will be made
available to the end users.  Assuming the end user have access to high speed
network and the parallel computer's cost effective computing power
is available to the general public, then terrain rendering is one of CPU
intensive applications that the user is interested.  To
capitalize on these features, a terrain rendering program for client-server
environment is needed.  However, the demand of computer graphics from the
consumer is high enough that mass production of specialize graphics hardware
is set in motion.  As a result, the use of parallel computer for
rendering graphics for general public is no longer cost effective.  Although the
custom program is no longer in favor, it gives us a in-depth
look on what the rendering techniques are needed.  This prove to be variable
knowledge for the design and implementation of VRML Terrain rendering.

In the following sections, we first provides an overview of the computer
graphics techniques and the details of our parallel rendering program.
Finally, we will presents our findings and improvements possible for our
program.

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

Computer graphics encompasses many techniques to turn a database of virtual
objects into images.  The objects can be defined analytically or
explicitly~\cite{Foley}. 
A collection of analytically defined elementary objects, such as, cone, sphere
and plane, is governed by the corresponding mathematical equations.
A complex analytical object can be formed by ``adding'' or ``subtracting''
using multiple elementary objects.
A explicitly defined object formed by a collection of 3D points (vertices)
and the interconnection of these vertices (edges).  Polygons are used to
emulate the surface, which is defined by the vertices and edges.

After the objects is defined, a wide range of rendering techniques can be used
to generate the final image.  Some of the techniques, such as ray tracing and
radiosity, are extremely computationally expensive, but they can be used to
generate photorealistic 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
acceleration hardwares.  However, the image generated is sometime ``appear''
artificial to the viewer.  This is due to the approximation used by the these
techniques.

\begin{figure}
\begin{center}
\begin{picture}(365,329)
\psfig{figure=rendering_pipeline.eps,width=365pt,height=329pt}
\end{picture}
\end{center}
\caption{A generic four stage rendering pipeline consists of transformation,
visible calculation, clipping and rasterization.
\label{rendering-pipeline}}
\end{figure}

We now presents a generic rendering process with a four stages rendering
pipeline for rendering explicitly defined 
objects (figure~\ref{rendering-pipeline}).  First, the
objects in the database (the virtual world) is transformed to the orientation
specified by the user. This orientation also known as the viewpoint of the
image.  The transformation process uses a collection of 
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 axis,

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

scaling
\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}
and
geometric projection $(M)$ matrices.
Although rotation and scaling of a 3D object can be done by using
$3 \times 3$ matrices, only $4 \times 4$ matrices can provide the translation
and the full range of geometric projections.
 
In order to use $4 \times 4$ transformation matrices, the coordinates,
$(x, y, z)$, of a 3D object is mapped into 4D coordinates, $(x, y, z, 1)$.
This maps the original 3D world into one of the 3D subspace
of a 4D universe.  The rotation and scaling is still operated
within the 3D subspace, so an identity matrix is used for the 4th dimension.
The operating principal for the $4 \times 4$
transformation process is to identifies the corresponding points in the other 3D
subspaces for all points in initial 3D subspace.  By projecting the points from
the other 3D subspaces back onto the initial 3D subspace,  special viewing
specification, such as, a perspective view, can be satisfied.  The perspective
transformation scales the objects inversely proportional to the distance, $d$,
between
the viewpoint and the objects. Hence, the far away objects in the virtual world 
appear smaller.

\begin{eqnarray}
M(d) & = &  \left[
\begin{array}{cccc}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 1 & 0 \\
0 & 0 & \frac{1}{d} & 0 \\
\end{array} \right]
\end{eqnarray}

An arbitrary view orientation requires a combination of the forth mentioned
$4 \times 4$ transformation matrices to put the objects in their proper
positions.
These matrices can be multiplied together to form one single $4 \times 4$
matrix
\begin{eqnarray}
U & = & T(x, y, z) R_{x}(\theta) R_{y}(\theta) R_{z}(\theta)
   S(s_{x}, s_{y}, s_{z}) M(d)
\end{eqnarray}
The combined matrix, $U$, is applied to each vertices of all the objects.  The 
transformation of one vertex requires 12 floating-point multiplication and 12 
floating-point addition.

Second, 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 reduce the total amount of effort needed to generate the
image.  A object is visible, if its volume intersect or contain within the
visible zone of the virtual world.  The intersection 
calculations can be CPU intensive, if the object have a complex shape.  A
transparent, simple
bounding volume, such as, a box or a sphere, which surrounds the object, can be
use to simplify the calculation.  Since the bounding volume is larger than the
object, some of the unseen objects will be marked as visible because their
bounding volumes are in the visible zone.  Tighter bounding volumes imply
a more accurate visible objects list.  A trade off is made between the
amount of work on intersection calculation and the amount of 
rendering effort saved, and the goal is to reduce the overall work to render the
image.

Third, all visible objects are projected onto a 2D plane that contain the
image area and the 2D objects are clipped according to the boundary of the
image area.  The process of projecting the 3D objects on to a 2D plane is
similar to the process of taking a photograph of the object.  For a explicitly
defined object, the projecting process consists of ignoring the $z$ coordinate
of all vertices and scale the 2D object to the size specified by the image size.
This also provides the first
opportunity of eliminating some the invisible polygons of the objects from
further 
processing.  Then the 2D projection is checked against the location and
size of the image area.  The invisible objects escaped from previous step
will be eliminated entirely.  New polygons are computed
by clipping the invisible portion of the partially visible
polygons.  The clipping computation is done by using
2D line equation and 2D line intersection calculation.

Finally, the 2D objects are drawn onto the image buffer pixel by pixel, which
is called rasterization.  At the pervious step, all the remained 3D objects are
projected onto the same 2D plane by ignoring the depth information, the 2D
objects may overlap each other on the 2D plane because the distance 3D
objects may be located behind the closer 3D objects.  An arbitrary order of
reaterizing the 2D objects may put the distance objects in front of a closer
object.  The correct order of rasterization can be retrieved form the $z$
coordinate of the 3D objects by a depth order search on the
list of visible objects.  Another popular method, which is called Z-buffer
method, is to record 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 color of the pixel is influenced by the properties of the light,
the properties of the object
and the approximations involved in computing the pixel color.
The properties of light consists of its position and color.  The relative
position of the light and the object determines the color intensity of the
pixel.  If the pixel is facing away the light, then, the color will be darker
(less intense).  The facing direction of the surface, where the pixel is
represented, is encoded by the surface normal.  A surface normal pointing
toward the light indicates that the surface is illuminated.  However, this does
not take into the account of shadows casted by the near by objects.  Omitting
shadow is one of the many approximations to minimize rendering resources
and time.  If the pixel is illuminated, then the color of the pixel is
a combination
of the color of the object and the color of the light, for example a red
object under a blue light appears black.

\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 shade of grey.
\label{Gouraud-shading}}
\end{figure}

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

The properties of the object not
only consists of its color and position, but also the surface properties of the
object.  However, rendering a realistic surface properties (all the groove and
bumps) is very computational
intensive, so approximations, such as Gouraud and Phong shading, are used to
compute the final contributing color of the object.
Gouraud shading is also known as intensity interpolation shading~\cite{Foley}.
The color of the object is first computed at the vertices of the object
Then the color of the pixel is computed by linear interpolation using
the near by vertices as shown in figure~\ref{Gouraud-shading}.
Phone shading is normal-vector interpolation shading~\cite{Foley}.
Instead of interpolating the color value of the vertices, the normals from the
near by vertices is calculated for the pixel (figure~\ref{Phone-shading}).
After the
normal is obtained, the color of the pixel is computed the same way as a
vertices.  Phone shading an object is obviously more computational expensive
than Gouraud shading because of additional work associated with computing the
normals.  However, Phone shading gives more realistically shaded object.
If the highlight, for example, falls between the vertices, then
the highlight will be missing from the Gouraud shaded object.
The highlight of a rotating Gouraud shaded object ``appears'' blinking
if the highlight alternately falls on and off vertices.
However, a Phone shaded object's the highlight will remain unchanged.

Further improvement can be extended from Phone shading algorithm by adding
slight perturbation to the computed normal.  This technique is called bump
mapping, which we have suggested using it in section~\ref{earth-bump-map}. 
We recall that the object is formed by a collection of polygons, which
represents a smooth surface.
The computed normal for Phone shading is obtained from the interpolation of
the vertices' normals of a polygon.  As a result, the normal varies uniformly
across the surface of the polygon.  However, if this normal
is perturbation before computing the highlight, then the color intensity will
varies accordingly.  Hence, the illusion of a bumpy is created.

% 1) texture map

\begin{figure}
\begin{center}
\begin{picture}(293,274)
\psfig{figure=scan-convert.eps,width=293pt,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 vertices list is indicated by the number.
Ths solid lines indicate the boundary of the screen pixel.
The dish lines are the scanlines. 
\label{scan-convert}}
\end{figure}

The rasterization of the 2D polygons is done by scan convert algorithm.
Assuming the polygon vertices are in clockwise order in the polygon's vertices
list (as shown in figure~\ref{scan-convert}).
The algorithm first find the starting scanline (scanline A), which is
determined by $\lceil y \rceil - 0.5$, where $y$ is the Y coordinate of
highest vertex in Y direction.
The left edge (edge defined by vertex 4 and 3) and the right
edge (vertex 4 and 1) of the polygon is locate on the starting scanline.
All the pixel on the scanline between the left edge and the right
edge is visited.  After visiting all these pixel, the scanline is moved
downward.
If the scanline landed on or just pasted a left or right vertex (scanline B and
C), another
vertex on the same side is found.  This process repeats until the left vertex
equal to the right vertex (scanline D).  Then all the pixel bounded by the
polygon is visited.

The image is completed after all the object is rasterized.
The entire process is illustrated
in figure~\ref{rendering-process}.  In the next section, we will showcase the
terrain rendering process and the computer graphics techniques use by it.

\section{Terrain rendering}

In this section, we will describe our terrain rendering process.  Although our 
terrain rendering process is modeled after the generic rendering pipeline,  we 
have optimized the process to capitalize on the data and viewing 
characteristics.

\begin{figure}
\begin{center}
\begin{picture}(365,329)
\psfig{figure=terrain-data.eps,width=365pt,height=329pt}
\end{picture}
\end{center}
\caption{The format of the DEM and Landsat images. \label{terrain-data}}
\end{figure}

The data for our terrain rendering program is composed from two datasets, the 
Digital Evaluation Model (DEM) and the Landsat sensor values.  Both datasets 
are coregistered to eliminate the visual distortion from the mismatch of the 
datasets.  The format of
the data is shown in figure~\ref{terrain-data}.  The header information of the
input files contain the UTM coordinates of the area and the size of data.  The
body of the data file is a 2 dimension array of integer values.  The Landsat
image is divided into 7 color bands and stored in 7 datafiles.  Six of the
seven bands is used to recreate the true color image.
The original DEM form USGS have a 100m resolution. Because the 
Landsat data have a 30m resolution, a 30m resolution DEM is constructed , so 
that there is a one to one mapping of the height value to the sensor values.
The 30m resolution DEM is created by replicating every data point in the 100m 
DEM three times as shown in figure~\ref{100to30DEM}.   The Landset sensor 
values passes through an ad hoc mapping function to form the RGB color 
value of the terrain data.   This mapping function is formulated based on the 
suggestion from the people at Jet Propulsion Laboratory (JPL), who helped 
process the Landsat data for the JPL terrain rendering program~\cite{}, and 
``looks good''' criterion.

In pervious chapter, we have apply wavelet to take advantage of the physical 
characteristics of the terrain data, so that the data in the visible and the 
neighboring region (the active region) is more manageable.   However, the 
generic rendering process call for  the transformation of the entire data set 
before the active region can be found by the visible calculation.  The 
transformation of a scene of 100km $\times$ 100Km (3333 $\times$ 3333 
vertices) terrain dataset requires 133M multiplication and 133M addition. A 
dataset, which covers the entire new York state, have approximately 10 scenes.  
Hence, the real time rendering by using the generic rendering pipeline of a 
large terrain dataset will be prohibitively expensive even for a large parallel 
computer.  As a result, we repackage the terrain data, so that the overall 
rendering work can be reduced.

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

The repackaging of the terrain data increases its modularity, so it can be
handle more flexibly.  The 2D terrain data is repackaged by dividing it into a
collection of 2D sub-blocks with the boundary overlapped as shown in
figure~\ref{repackage}.  This conformed to the suggestion in
section~\ref{model-joining} on joining the terrain models.
Finally, the data is stored in a vertices list, an
edge list and a polygon list. Each entry of the vertices list contains the 3D
coordinates and the RGB color of the vertex.  The edge list entry have
two indices of the starting and ending vertices.  The polygon list entry
consists of three indices of its edges.  This storage format is show in
figure~\ref{entry-storage}.
After the repackaging, the original terrain data set is broken into many
computer graphics
objects.  This allows us to manipulate the terrain data set in the $O(M)$,
where $M$ is the number of objects, instead of $O(n^{2})$, where $n$ is the
size of the terrain data.  Although these objects are connected, they are
independent from each other because there is no edge crossing the package
boundary.  As a result, elimination of an object will not introduce
inconsistency, such as incomplete polygon or undefined vertices.

In addition to the repackaging of the data, the viewing characteristic allows 
shading to be done in the preprocessing step because of the invariance of the 
dataset and the best viewing light condition.   The terrain features are seldom 
change 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 
survey for the maps are done at the beginning of the 20th century.  Because of 
the satellite orbit and natural obscures, such as, cloud cover and fog, a clear 
satellite image is a chance encounter.  Hence, the best Landsat images is 
collected over the years, and the image for our test area is 5 years old.  This 
foster an 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

The pre-select sun position to provide the best viewing condition also enable 
the shading of the terrain to be computed in the preprocessing phase. We 
assume the user is in a dry and clear summer day with 12 hours day light, and 
the sun swaps across the sky form east to west in a $180^{\circ}$ arc.  The sun 
can only provide the best viewing light condition in a limited range of angle.  
For the 2 hours period around dawn and dust ($30^{\circ}$ just above the 
ground), the shadow of the terrain is elongated. The shading techniques, which 
is used for rendering the terrain data, does not provide a good approximation 
for the elongated shadow (see section~\ref{} for more detail discussion).  For 
the 2 hours around noon ($\pm 15^{\circ}$ of the vertical), there is not enough 
shadow to provide a quick navigational reference.  For example, without any 
shadow, a hill on the east of the viewer looks exactly like the hill to the
west.  Hence, the user is easily lost within the virtual world.  These three
time periods (total of $90^{\circ}$ arc interval) are not suitable for viewing
the terrain data.  The rest of  time the sun is either angled slight to the east
or to the west (two 
$45^{\circ}$ arc period as show in figure~\ref{global-angle}).  The rendered 
images under both lighting angle periods are similar except that the shadow is 
pointing at a different direction.  Therefore, the global light condition can
be pre-selected to provide the best viewing condition, and the shading does not 
need to be done at runtime.  Since the preprocessed terrain data is being
reused by the rendering program, the repackage operation and the shading needs
only to be performed once.   Consequently, an independent preprocessing
pipeline is added to the terrain rendering process.

The entire terrain rendering process not only gain an addition pipeline, but
the visible calculation is also splited into two stage to exploit the advantage
offered by the repackaged data (figure~\ref{terrain-pipeline}).  In
section~\ref{visible-calculation}, the coarse visible calculation is describe
as rasterizing the 2D projection of the view volume (bounding polygon of all
visible data) onto the terrain data grid..  All the grid sections that are
with in or partially within the bounding polygon are assume to be visible.
Since the viewer orientation is limited (see section~\ref{visible-calculation}),
the 2D view volume projection onto the data plane is always triangular.
Because of the limited size of terrain data grid the bounding edges of the
grid promotes the triangular projection to a higher edge count polygon.  There
are total of three possible polygon shape without considering the rotation
symmetry as shown in figure~\ref{view-projections}.  The limited number of
configuration allows us to hardwire the rasterization to speed up the
computation.  The computation consists of maximum of $2n + 4$ number of
intersection calculation.  


We will discuss the detail of  our parallel implementation of this rendering 
pipeline.

\section{Parallel terrain rendering}
\label{parallel-rendering-pipeline}

   In previous section, we have provided an overall view of the terrain
rendering pipeline.  We now discuss the design and implementation
of a parallel terrain rendering process.  The motivation of a parallel
implementation of the terrain rendering is to handle the large amount of data and
to increase rendering speed without using specialized hardware.

When the DEM is repackaged into computer graphics objects, the total
storage space needed increases by a
factor of 15.  The increase comes form the explicit definition of the height value
in 3D world, the encoding of the edges and polygons.
Although not all the data is needed at the same time, the data is still load in core.
As a result, a parallel computer offers enough aggregate memory to house the
data.

The large data size also indicates that there are a lot of polygons to render.  Since
each terrain grid section is independent form each other, the rendering can be speed
up by using multiple rendering pipeline line.

In order to increase the rendering speed and the size
of the data sets that can be handled, the terrain rendering program 
must be implemented in parallel.  The parallelism of the terrain rendering is
in the from of object parallelism, pixel parallelism and frame parallelism.
The object parallelism derives from the terrain data's characteristics of
non-overlapping and having
an unambiguous depth order.  Consequently, polygons can be rendered in
any arbitrary order and the characteristics ensure the same image is
generated.  As in all other graphics applications, all pixel in an image
are totally independent from each other.  This implies the pixel can also be
rendered in arbitrary order.
By combining object and pixel parallelism characteristics, the terrain
rendering program can distribute and render the data efficiently in parallel
without compromising the image quality.
Because the terrain program is interactive, the program has to generate
multiple images per second.  However, each image is generated with no
dependency on the pixel values of the previous images.
Consequently, the terrain rendering program can render a number of images
for different time steps in parallel.

Since the preprocessing is performed on the entire database, it is an good 
candidate for parallelization
\subsection{Parallelism}
\subsection{Data Structures}
\subsection{Preprocess}
\subsection{Runtime}

\section{Findings}
The goals of implementing a parallel terrain rendering program is to
increase the rendering speed and the size of the data that can be handled.
Rendering speed is quantify by the number of frame per second and the number of
polygon per second.

\section{Improvement}




