Glenn Elliott
Brian Manson
Remigio Parayao
12/02/04
CS 580 - Terrain Polygon Reduction
Introduction
Terrain data is abundant and in high detail from survey projects performed both through traditional and space-based methods. With so much survey data for relatively small areas, rendering this data in visualization systems can be quite a task and often overpowering computationally through brute force methods. To overcome being overwhelmed with data, elegant solutions attempt to reduce or generalize data in terrain areas where high details are not required. There are many ways to approach this problem and our technique is one of many developed in the field. Our approach is to generalize all the terrain data with spline equations, which can then be re-sampled at adaptive rates.
This report is divided up into four sections. The first discusses our source of terrain elevation data. The second section covers our technique of using quad-trees to subdivide the large terrain dataset into patches of varying size and detail. The third section discusses the math behind spline calculation. The last section covers implementation details and the difficulties we faced during implementation.
Digital Terrain Elevation Data (DTED) Background
In support of military applications, the National Imagery and Mapping Agency (NIMA) has developed a standard digital dataset format called Digital Terrain Elevation Data (DTED) which contains a uniform matrix of terrain elevation values. This data provides basic quantitative information for systems and applications that require terrain elevation, slope, and/or surface roughness information. We use Digital Terrain Elevation Data (DTED) as our source of terrain data in this project.
DTED data comes in various forms and are classified into five levels of detail. DTED Level 0 data is sampled at 30 arc seconds with a resolution of nominally one kilometer per point. DTED0 is derived from DTED1 and freely available on the Internet. DTED Level 1 has a resolution of 3 arc seconds, roughly 100 meters. DTED Levels 2, 3, 4 and 5 have resolutions of roughly 30, 10, 3, and 1 meter respectively. DTED 1 and 2 is available for purchase on the Internet. Levels 3, 4, and 5 are government controlled.
Our project uses DTED1 data (we found that DTED2 involved too much spline calculation time than we were willing to spend). The data extraction process starts by reading the header of the DTED file. The DTED file has information on the number of latitudes and the number of longitudes, in other words the header describes the resolution of the terrain map. Once the terrain size has been determined, the terrain elevation (which is in meters) can be extracted by extrapolating the latitude and longitude starting at the lower left corner of the grid. For every elevation points, we need to figure out the x and the y values. Our data modeler uses the extracted DTED, convert the latitude and longitude in meters using the Soldano inverse method. This method takes into account the curve and irregularity of the earth. The final data is a matrix of size (number of longitude * number of latitude). For each matrix cell it will have an x, y and elevation in meters.
Terrain Analysis and Quadtree Utilization
Fundamental to our approach in terrain polygon reduction is the ability to generalize large simple terrain patches like smooth slopes and plains while retaining the required detail of rough mountainous regions. Our approach to this was to recursively subdivide the terrain dataset into patches until each patch was of a reasonable smoothness. This was done through the use of a quad tree and our algorithm is as follows: Starting with the entire terrain dataset, a measure of terrain roughness is made. If this roughness is greater than an acceptable limit, then the patch is cut evenly into quarters and the process is repeated on each subdivision. This continues until all leaves in the quad tree data structure are of an acceptable smoothness or have reached a minimum size. The result is a quad tree structure of varying depths. Leaves containing large smooth areas of the terrain are higher up in the tree and then highest detailed areas are deep in the tree.
Generation of the quad tree is simple and straightforward given one has a method to measure terrain smoothness. We found three different methods for analyzing and classifying terrain smoothness and are as follows.
Standard Deviation of a Uniform Distribution –The premise behind this method is that areas that are simple and smooth are flat and areas that are rough, are not. To measure this we can calculate the standard deviation of the terrain elevations of a given patch from the uniform distribution. The equation for this is simple:
However, this approach is naïve. The standard deviation does not consider the locality of each elevation to those around it. As a consequence, is that a large smooth simple slope still has a high standard deviation around its mean. A suitable method must consider this case.
Standard Deviation of a Uniform Distribution with Fitted Plane – This method takes into account sloping plains by estimating the general slope of a terrain patch and then subtracting this sloping component from the patch dataset. This effectively flattens the terrain and allows us to calculate the standard deviation of the data from the uniform distribution. This approach is conceptually straightforward, however, the estimation of the fitted plane is difficult and requires quite a bit of computation.
Figure 1: The above depicts the 2D concepts behind the fitted plane method. A fitted line is estimated from the rough terrain. This line is then subtracted from the source. The sloping component is now removed and standard deviation calculations can be made.
Roughness Index – The roughness index method is an alternative to the fitted plane method in that it also accounts for terrain slopes. The roughness index is derived from an absolute measure of the second derivative of the terrain data. Simply put, it is a sum of the total (absolute) curvature of the terrain data. The equation for this is as follows:
Because the equation measures curvature, slope does not need to be considered. We measure the change in slope from point to point across the terrain patch not slope itself. While this method is simpler then the fitted plane method, it does have one great draw back. It requires that the terrain dataset be organized into a regular grid. This method cannot be used on scattered datasets while the fitted plane method can.
In the end, we implemented roughness index method. This is because our source DTED data is organized into a regular grid. If we were to extend our project to support scattered data, the fitted plane method would have to have been used.
Below are the results of our quad tree generation. Note that large patches cover the smooth regions while the mountainous areas are represented with smaller patches.
Figure 2: Small patches cover rough mountainous regions.
Figure 3: Large patches cover the smooth valleys at the mountain base.
Spline-Based Interpolation
In generalizing terrain data organized in a uniform grid-- two problems exist. First, how does one interpolate between points in the terrain? The result must manufacture data that has been excluded from the dataset. Furthermore, the interpolation must mimic the nature of the data (rocky, but sometimes smooth terrain) and the interpolation at the edges must blend well with neighboring points. Secondly, after the interpolation is complete, how can we find points that don't exist on any of the data axes?
Spline interpolation solves the aforementioned problems well for many reasons. The nature of the data is somewhat inconsistent, but often the points in the dataset belong to sloping curved surfaces (see figure 5.) Linear interpolation of points would create a result similar to the picture in figure 1, which would result in a very spiky and unrealistic terrain surface.
Figure 4: linear interpolation Figure 5: spline interpolation
However, the decision to solve the interpolation problem using splines is not sufficient- we must sort through the myriads of spline classifications that the math world has to offer. Ordinarily, splines are classified by the smoothness of their continuity. Continuity is defined by the similarity of the derivative at the joining point of two curves. Figure 6 shows an example of C0 continuity where it is sufficient to simply share a point at the joining spot of two curves.
Figure 6: C0 continuity Figure 7: C1 continuity
Figure 7 shows an example of spline curves that shares the same first derivative at the joining point. Note that the two curves need not join at the location of the middle point. Point sharing is not a requirement of C1 and above continuities. However, some splines do allow joining at the middle point, as well as, sharing a higher order derivative. C2 continuity takes smoothness to next level by sharing second order derivatives at the joining point (see figure 8.)
Figure 8: C2 continuity
As mentioned earlier, our continuity wish list includes two facts: point sharing at the middle point so that every location in the dataset is properly represented, and a higher order derivative similarity at the joining point to ensure a maximum smoothness. Fortunately, there is a spline that satisfies both of these criteria- the Catmull-Rom spline.
Ed Catmull, the co-founder of Pixar, and Raphael Rom, developed the Catmull-Rom spline in 1974. The Catmull-Rom spline, or CR spline as it is sometimes called, was an improvement on other splines of higher order continuity because it guarantees that all the control points will pass through the curve. The most popular predecessor to the CR spline is the Bezier spline, which holds a C2 continuity, but cannot guarantee point sharing at the middle point. The problem that Bezier splines create is that with dozens of points per spline the programmer must choose between alternating C0 and C2 continuity by joining multiple splines, or consistent C2 continuity with complex methods of inter-spline blending.
The Catmull-Rom spline is simple to compute, however, it introduces some new problems into the terrain interpolation procedure. The computation of a CR spline involves a minimum of 4 control points and a t that determines an offset between curves (see figure 6.) The nice thing about CR splines turned out to be that neighboring spline equations fit together seamlessly because of the use of future points and previous points.
The Catmull-Rom computation is as follows:
Note the role of t in the equation above and how it relates to the curve offset in the spline from figure 9.
Figure 9: Catmull-Rom spline
Notice that as we loop through a series of points, we must observe the second point as the current point. The number of points mean that we must specify a point prior to the current one and two points beyond the current point; herein lies a substantial difficulty. Often, we will be on one edge or another of a patch that needs to be interpolated. My first solution was to take any two points that are given, and generate a future and/or previous point by computing the slope between the two points and estimating the probable location of another point. More often than not, the guess was incorrect, and adjacent patches would be disjointed. The best result was given by reading the elevation data outside the current patch to serve as future/previous points. The only downside to this approach was that we would have to disregard a 2 point perimeter of the dataset. Given the fact that our dataset had hundreds of thousands of points, this was a small sacrifice.
The final problem was the determination of intra-spline points after the creation of the spline network (see point p in figure 7.)
Figure 7: spline patch Figure 8: bilinear interpolation result
There are many methods of retrieving point p from a patch of splines- the most obvious being bilinear interpolation. However, this did not create the most desirable results (see figure 8.) Figure 8 illustrates how bilinear interpolation has ruined the integrity of one patch. The CR spline representation of points lying on a line is mostly flat with some minor curvature. The bilinear representation of this surface accentuates the minor curvature of the surrounding splines. The best result came from taking a corresponding point on each surrounding edge and averaging those four points. The result was much smoother and more natural.
Project Implantation
The following diagram shows the flow of data in our project from beginning to end.
Figure 9. The above figure depicts our project data flow. The data modeler will extract the DTED and create a matrix of size specified by DTED. Each matrix cell will have x, y and elevation. The final matrix will then be feed to the rendering function and ultimately produce a 3d terrain.
The process was broken up among the three of us into DTED extraction, quad tree generation, spline calculation, texturing and rendering. DTED, quad tree, and spline topics have already been discussed. For rendering, we decided to use the DirectX 9 graphics API. This allowed for us to concentrate on the calculations required to support our project aims and not worry about the details of polygon rendering. Further, we also benefited from hardware accelerated graphics computation. In our implementation we often found ourselves working with more than 5 million triangles at a time. Without hardware acceleration, we would not have been able to handle so many polygons.
Rendering was simple to do. Each quad tree leaf sampled the spline equations for its patch at an adaptive rate determined by the leaf’s depth in the quad tree. The set of points were then used to form a mesh. Normal calculations and texturing were then performed on the mesh vertices. This data was then pushed to DirectX as an indexed triangle list to be rendered.
Project Difficulties and Lessons Learned
Difficulties
The project began with a DTED dataset and the potential to show some photo-real results given the intricacies of the data. After our first group meeting, we had come to the conclusion that simply representing the data in a DirectX application was not enough to constitute "research." For the next week or two, we struggled to find a direction of study using the DTED data as our starting point. I think it's a fair assumption that most research doesn't start with the data. Nevertheless, it has been an interesting experience. After several Google sessions, we found a very complicated paper regarding mesh smoothing using quad trees and Bezier splines. We understood some of the broad points, but we were generally confused about the content of this paper. After some discussion, we came to the conclusion that (given the few things we understood about mesh smoothing) we could attempt this topic as our project. We diverged from many of the key points of the paper, and placed our emphasis on polygon reduction rather than smoothing.
Our main bug in the process of building the project related to DirectX rendering. For several hours, we combed the code for an explanation of why not all the polygons in each patch were being rendered. The correct number of polygons seemed to be specified and allocated. We stepped through the creation process of the patch mesh and could find no errors. If the terrain wasn't being sub-sampled or super-sampled everything worked fine. After a long period of debugging, we found that the renderer was being told to render only the unadjusted amount of polygons. This explained why in sub-sampling some polygons would render to create overlapping patches, and in super-sampling polygons would be missing altogether. Having debugged a DirectX application now, I can say with certainty that it is not an API that is conducive to easy debugging. The majority of the code is very low-level and hard to trace.
Lessons Learned
Making a demo-able product can be an art form. We feel our final product fit more into the category of useful and innovative, than easily presentable. The worst part of the executable is the enormous amount of time that it takes to start. There is a lot of precomputation that takes place upon startup. After the initial startup, the time involved in flying through the scene could be titanic depending on the rate of sampling. For the most part, frame-rates were above expectation. It would have been worthwhile to create some interface for increasing/decreasing the sample rate during runtime so that people could see the transformation, or lack thereof, in the terrain. This feature was on our wish list, but we were not able to reach it in time.
Another feature, which would have improved our final product, was the use of the patch statistics in spline computation. The patches were formulated with a generic spline algorithm. Given a grade of "flatness" in the patch, the spline could have been blended with a plane or calculated with some other approach to spline interpolation. Given a patch that correctly accounts for flatness, bilinear interpolation could have been used to generate more accurate surfaces.
Lastly, we were pleased by our ability to segment the project into 3 even parts, but our planning lacked some specifics that we found out about later in the coding phase. If we had discussed the inputs/outputs of some of our functions, as a group, we would have benefited from not having to re-write small parts of functions. There was some misunderstanding at one point regarding whether the spline code should access areas of the terrain by an i/j index or the latitude/longitude. The amount of code re-written was minimal, but it did slow down our progress a bit.