Subsurface Characterization with Three-Dimensional Natural Neighbor Interpolation

REFERENCES

Abstract

While recognizing the merits of many of the numerous trivariate interpolation schemes [Alfeld,84] [Alfeld,89] [Davis,86] [Dahman,89] [Dyn,90] [Franke,82] [Franke,80] [Lancaster,86] [Petrie,87] [Rescorla,87] [Watson,84] [Watson,85], most do not provide the combination of features needed for effective interpolation for site characterization. One method, satisfying these requirements is known as Natural Neighbor Interpolation. An implementation of natural neighbor interpolation in both two and three dimensions is the subject of this paper. This scheme, first introduced by [Sibson,81], is an inherently local scheme that performs very well with clustered data [Watson,87]. While Sibson addressed interpolation only, a modification of natural neighbor interpolation is also presented in this paper which allows for extrapolation.

Although some literature is currently available on the subject of natural neighbor interpolation [Sibson,81] [Watson,87], very little is presented about the specifics of implementation and no previous work, that we are aware of, directly addresses implementation in three dimensions. The objective of this research is to define a trivariate algorithm for natural neighbor interpolation that incorporates some of the unique requirements of the required application. For most interpolation schemes, extending their applicability from two dimensions to three, is a trivial matter. Since Natural neighbor interpolation involves Dirichlet and Delaunay tesselations, its extension to three dimensions is not as trivial. For this reason, a two dimensional prototype is first developed and then extended to three dimensions. We first review the theory behind natural neighbor interpolation followed by a discussion of the details of implementation for both two and three dimensional data. Extrapolation using the methods defined for natural neighbor interpolation is also discussed as well as a limited number of sample cases tested using an implementation of the scheme.

2. Natural Neighbor Interpolation

Figure 1. The Delaunay and Dirichlet tesselations of a set of scattered data

Natural neighbor local coordinates are illustrated in figure 2. Points 1-10 are scattered data points where the z value is known and P_n is a point where we would like to interpolate a z value. Temporarily inserting P_n into the Delaunay tesselation will cause the Dirichlet tesselation to change resulting in new Thiesson areas for the polygons in the neighborhood of P_n. The dashed lines show the Dirichlet tesselation before P_n is temporarily added to the TIN and the solid lines show the Dirichlet tesselation after P_n is added. Only those data points whose Thiesson polygon are changed will be used to compute the interpolant. In this case only points 1,4,5,6 and 9 will be used. The local coordinate, l_m(n) for each of these points with respect to P_n is defined as the area shared by the Thiesson polygon defined by P_n and the Thiesson polygon defined by the respective data points before P_n is added. The greater this common area, the larger the resulting local coordinate and the greater the influence on the interpolant.

Figure 2. The Dirichlet tesselation of a set of scattered data before and after P_n is inserted.

It is apparent that l_m(n) will be a number between zero and unity. If P_n is at the same location as P_m then l_m(n) will be unity. The local coordinate also has the property that the sum of all local coordinates l_m for any interpolant P_n will be 1.0. Similar to inverse distance weighted schemes, in general, the greater the relative distance P_m is from P_n, the smaller its influence on the final interpolated function value. Unlike inverse distance weighted interpolation, however, the weights are also dependant on the spacial relationships of the influencing data points with respect to each other.

The neighboring data points involved in the interpolation of P_n are those points defining triangles in the Delaunay triangulation that are adjacent to the temporarily inserted P_n. If we define k(n) as the Thiesson area of P_n and km(n) as the difference in Thiesson area of neighboring data point, P_m, before and after P_n is inserted, then the local coordinate lm(n) is defined as:

(1)

Knowing the local coordinates or influence of the neighboring data points to P_n, the function value may be computed by summing the function values at neighboring data points weighted by their respective local coordinate as follows:

(2)

where k is the number of "natural neighbors" to the interpolant P_n.

Additional refinements may be made to this equation to incorporate gradients and blending functions. Sibson suggests the following formulation for the interpolant:

(3)

where b_m is the estimated gradient at P_m and u_m and u_n are the coordinates at P_m and P_n respectively. The blending function, B, is a function of the local coordinate, lm. . The equation used for the blending function in this case is a simple hermite polynomial blending function:

(4)

The blending function has the effect of enforcing C² or second derivative continuity on the surface.

Having defined the natural neighbor interpolant, there are several issues that must be resolved in order to implement such a scheme. These are divided into the following categories:

1. Defining the Delaunay tesselation. In order to define the local coordinates for a data point, its Voronoi polytope must be defined. The preliminary step to this is first to define the Delaunay tesselation for all data points. The compliment Dirichlet tesselation comprising the network of Voronoi polytopes can then be defined
2. Defining the Thiesson polytopes and their areas or volumes. The geometry of each Voronoi polytope is defined by using the circumcircle or circumsphere centers of simplices comprising the Delaunay tesselation. Algorithms for computing the area or volume of the resulting polytope are presented. Local coordinates, lm(n) , for each neighboring data point can then be computed from equation (1).
3. Estimating gradients. Methods for gradient estimation are presented that are consistent with Natural Neighbor interpolation that will improve the continuity of the resulting function.

Once each of these components has been assembled, equation (3) can be used to compute the final interpolant. Individually, each of these areas could be the topic for a lengthy discussion. Rather, we will only touch briefly on some of the specifics of implementation as they pertain to natural neighbor interpolation and its application to subsurface characterization.

3. Defining the Delaunay tesselation

In order to define the neighborhood and local coordinates for each data point in a set, natural neighbor interpolation requires that the data points first be tessellated into a network of simplices that maintain Delaunay criterion. It also requires that the simplices define a convex hull. The Delaunay tesselation of a set of data points provides the basis for defining the Dirichlet tesselation and the resulting area or volume of influence for each data point.

Methods traditionally used for Delaunay tesselation [Baker,89] [Bowyer,81] [Lawson,77] [Lee,80] [Watson,81] involve first defining a bounding simplex or convex hull of existing data points. Simplices are subdivided as successive data points are added and edges or faces are checked and swapped if necessary to maintain Delaunay criterion. The Delaunay criterion requires that vertices of neighboring simplices not be contained within the circumcircle or circumsphere of individual simplices. These methods typically require that the entire data set be defined before beginning the tesselation process. As each data point is added, it is added within the existing convex hull of simplices. In addition, some methods do not guarantee a convex hull once tesselation is complete. Because of the unique requirements of the application, a method for Delaunay tesselation was defined that would allow points to be added either inside or outside of the existing convex hull as well as maintaining the convex hull during the process of point insertion. This provided the flexibility of incrementally adding data points without prior knowledge of the boundaries of the entire data set.

Within the context of this application it was found that the tesselation procedure could be used for two different purposes. First, the data points are tessellated to maintain Delaunay criterion, providing a basis for defining the Thiesson areas or volumes. Secondly, in order to compute the relative weights for neighboring data points, a ratio of Thiesson areas before and after an interpolation point is temporarily added to the network is computed. Temporarily adding the interpolation point, also requires retesselation. Using an incremental approach to tesselation is also useful in the latter context.

The first step in the tesselation of a set of scattered data points to maintain the Delaunay criterion is to define an initial simplex. The initial simplex is defined by selecting vertices from the known data points. Ideally, the data points should be as widely separated as possible, since the point insertion process is faster for interior points than for exterior points. The initial vertices should be non-coplanar for the 3D case and non-collinear for the 2D case. After defining an initial simplex, each subsequent data point is individually inserted into the Delaunay tesselation. By locating each point with respect to the current tesselation, new simplices are added within the data structure as the tesselation grows.

One of the most important features of this tesselation algorithm is that the perimeter of the Delaunay tesselation will always be the convex hull of the data points. This characteristic is essential when defining Thiesson polytopes. It guarantees that any Thiesson polytope in the Dirichlet tesselation is convex and ensures that boundary Thiesson polytopes will be open-ended. Taking advantage of these properties is important to the algorithms that follow.

3.1 Point Location

3.1.1 Point Insertion: Adding vertices inside the convex hull

When a new insertion point is found to be contained within the convex hull of existing points, a technique described by [Lee,80] is used to insert the point. This method involves finding all neighbor simplices to the new insertion vertex (figure 3a). Neighboring simplices are defined as those simplices whose circumcircle or circumsphere contain the point. This list is compiled by propagating outwards from the triangle containing the point and recursively checking adjacent simplices. A list of boundary edges or faces is then compiled which describe the perimeter of the neighboring simplices. Finally the neighboring simplices are deleted (figure 3b) from the network and the boundary edges or faces are connected to the new insertion point to define the new tesselation (figure 3c).

	Figure 3a. Initial TIN showing where the new point is to be inserted. A list of neighbor triangles to the insertion point is compiled.
	Figure 3b. A boundary edge list of the neighbor triangles is compiled. The neighbor edge list is then deleted leaving a void in the TIN.
	Figure 3c. The void is now retesselated by connecting the boundary edges to the new insertion point

Figure 3. Point insertion inside the convex hull

3.1.2 Point Insertion: Adding vertices outside the convex hull

Because Natural neighbor interpolation requires that the boundary of the TIN be the convex hull of all of the data points, any point inserted on the outside of the TIN must form simplices with all points on the boundary of the TIN that are visible to it. An edge or face on the boundary of the TIN is defined as visible if, for example, an observer standing at the location of the insertion point has a direct line of sight to the edge or face without being blocked by any other simplex in the network. Mathematically, this is defined by determining the perpendicular distance from the point to the infinite line or plane of the boundary edges or faces. If simplex boundaries are always defined with inward pointing normals, any negative distance will define a visible edge or face. The edges that would be visible to a new insertion point are shown in figure 4a.

	Figure 4a. New point is to be inserted into the TIN. All visible edges to the new point are found and a boundary edge list is compiled
	Figure 4b. Any adjacent triangles to the boundary edges whose circumcircle contains the new point are deleted from the TIN. New visible edges are added to the edge list.
	Figure 4c. New triangles are created by connecting each edge in the edge list with the new point.

Figure 4. Point insertion outside the convex hull

Similar to the algorithm to add a point inside the TIN, for the two dimensional case a boundary edge list is compiled. This list will contain all edges that are visible to the new point. New triangles can be formed by connecting the new insertion point with each boundary edge in the list. The boundary edge list, however is not sufficient to generate new triangles that maintain the Delaunay criterion. All adjacent triangles to the edges in the boundary edge list must be checked to ensure that their circumcircle does not contain the new insertion point. If the circumcircle of a triangle is found to contain the new point, it is deleted from the TIN and its two opposite edges are placed on the boundary edge list. These two new edges in the list are then checked to see that their adjacent triangle does not contain the new insertion point. Figure 4b illustrates the procedure for checking boundary edges.

Once the final boundary edge list has been compiled, the new triangles may be generated. New triangles are formed by connecting the new insertion point with each edge in the edge list. The final tesselation with the newly inserted vertex is shown in figure 4c.

While several algorithms for Delaunay tesselation in three dimensions have been presented [Baker,89], the method described here seems to be most conducive to the application and is an extension of the two dimensional case. Once it has been determined that the new insertion point is on the outside of the existing TIN, like the two dimensional point insertion algorithm, a list of perimeter tetrahedron faces visible to the point is compiled.

In order to create the desired list of visible faces, we start from a tetrahedron face visible to the exterior point and recursively check adjacent boundary faces until all visible faces have been found. Since the TIN, throughout the tesselation process, is guaranteed to maintain the convex hull of the data points, no concave regions will exist. This is significant because it means that once a face has been determined to be not visible to the point, any further recursive checking of adjacent perimeter faces is unnecessary.

Once a list of visible faces is compiled, the tetrahedron containing the faces must be checked to ensure that their circumspheres do not contain the insertion point. Upon determining that a perimeter tetrahedron contains the point, the tetrahedron is deleted from the TIN and its perimeter face is removed from the perimeter face list. By deleting the tetrahedron, new faces will be exposed. These exposed faces must now be added to the perimeter face list and, in turn, checked to see that their adjacent tetrahedron circumspheres do not contain the insertion point.

The final step in the three dimensional exterior point insertion algorithm is to include the new insertion point in the tesselation. New tetrahedra are created by using each face from the perimeter face list and the exterior insertion point to define the tetrahedra.

The point insertion procedures are also used for an additional purpose. The determination of the local coordinates for a particular interpolation point requires that a ratio between the areas or volumes of Voronoi polytopes before and after a point is temporarily inserted into the Delaunay tesselation be computed. While this procedure can also be used for temporarily inserting a point, since the point is only briefly inserted into the network for puposes of local coordinate computation, the simplices deleted from the TIN cannot be permanently removed from the network. This implies that the removed simplices be saved and replaced after the local coordinates have been computed.

4. Defining Thiesson Polytopes

As previously mentioned, Natural Neighbor interpolation requires that the area or volume of Thiesson polytopes (also known as Voronoi polytopes) be defined. Figure 1 shows that for the points at the perimeter of the TIN, their Thiesson polygons have infinite area. For this reason a boundary is defined to limit the size of the polygon. In practice, the boundary would be a property line, boundary of a drainage basin or any other convenient polygon that will enclose all of the data points. For this study we will limit the boundary polygon to a simple rectangle aligned with the coordinate axis. This will simplify some of the calculations, however with small modifications, the algorithm may be generalized to use an arbitrary boundary polygon provided it remains convex. Likewise, a right parallelepiped will be used for the boundary in the three dimensional case.

Once all points have been tessellated to maintain Delaunay criterion, it is a simple matter to define the Dirichlet tesselation. [Bowyer,81] shows that for a data point in a TIN, the vertices of its Thiesson polygons are defined by the circumcircle centers of its adjacent triangles. In addition, the edges of its Thiesson polygon are the perpendicular bisectors of the Delaunay triangle edges. Similarly, for the three dimensional case, it can be shown that if a TIN has been tessellated to maintain the Delaunay criterion, the circumsphere centers of adjacent tetrahedron to a data point will define the vertices of its Thiesson polyhedron. The faces for the Thiesson polyhedron are defined as the perpendicularly bisecting planes for all tetrahedra edges radiating from the defining data point.

There are two cases to consider when determining the area or volume of a Thiesson polytope. The first and simplest case is when the entire polytope is contained within the window defined by the bounding rectangle or parallelepiped. This occurs when the defining data point is not on the perimeter of the TIN. Voronoi polytopes for this case are always finite. The second case occurs when the defining data point is on the boundary of the TIN. In this instance the polytope is open-ended or infinite and must be clipped by the window. One additional case arrises when a portion of a finite Voronoi polytope falls outside the bounding window. This occurs when the circumcircle or circumsphere centroid of an adjacent simplex to the defining data point falls beyond the boudaries of the window. Since the polytope must be clipped, this case can be treated along with the latter case. The two situations will be treated separately.

4.1 Polytopes Contained Within the Window

The vertices for a Thiesson polygon on the inside of the TIN are defined by the circumcircle centers of adjacent triangles. For example, the area of a Thiesson polygon defined in figure 5 is defined as follows:

(5)

where n is the number of vertices on the Thiesson polygon.

Figure 5. Thiesson polygon defined on the inside of the TIN

For the three dimensional case we take advantage of Stokes and Gauss's theorem [Goldman,91] to compute the volume of the Thiesson polytope. Stokes theorm is first used to compute the area of each face of the polyhedron and then combined using Gauss's theorm to comute the final volume. Stoke's theorm is stated as follows:

(6)

where Pk are the vertices of the polygon and N is the unit vector perpendicular to the plane where the polygon lies. For a polyhedron with planar polygonal faces S₀,...,S_n, Q_j is any point on S_j and N_j is a unit vector normal to S_j, then the volume of the polyhedron is:

(7)

Defining the polyhedron faces, S_j requires that each tetrahedron edge radiating from the data point be processed. The perpendicularly bisecting plane for each radiating edge will define the polyhedron faces. Vertices of each face, Pk are found by pivoting around the radiating edge. The circumsphere centers of tetrahedron sharing the common radiating edge define the vertices of a planar polygon. Finally, the normal N_j required for equation (6) is, by definition, the normalized vector defined by the outward pointing radiating edge. Combining this information, (6) may be used to compute the area of a face of a Thiesson polyhedron.

If we define Q_j as any of the polygon vertices of S_j, then knowing Area(S_j) and N_j, it is a simple matter to compute the final volume of the Thiesson polyhedron using equation (7).

Figure 6. Thiesson polygon defined on the outside of the TIN.

If, for example we want to find the area of the Thiesson polygon shown in figure 6, each of its candidate edges is determined and stored. Figure 7 shows all of the candidate edges for this polygon. The equations for the candidate edges are determined by looping around the data point in a counter-clockwise direction with the parametric direction of the line increasing to the right. In addition, the lines representing the boundary or window of the interpolation area are added to the list. These are inserted in order to clip the otherwise open-ended polygon. For an arbitrary shaped, convex bounding window, the equations for each edge of the domain would be computed and added as a candidate edge.

Figure 7. Candidate edges for a Thiesson polygon.

The next step in this procedure is to determine for each candidate edge, the intersection point with every other edge in the list. The intersection points are the candidate vertices of the Thiesson polygon. Parallel candidate edges are ignored.

Figure 8. Table to keep track of inside/outside information for candidate Thiesson polygon vertices.
(I=inside, O=outside, P=parallel)

Once the intersection points have been determined, a pass is made through each intersection point to determine whether it is “inside” or “on” all of the candidate edges. Figure 8 shows the table that is used to keep track of the inside/outside information. The sign of the distance from an intersection point to every edge is determined. If the distance from an edge is negative, it is marked as outside. After a complete pass through each intersection point, the remaining points are the vertices of the Thiesson polygon. The table shows the results of this procedure. “I” indicates that the point is inside or on all lines, “O” indicates that the point is outside at least one of the lines, and “P” indicates that the two lines are parallel.

To compute the final area of the Thiesson polygon a pass is made through the candidate edges. The area under each edge is determined between two non-equivalent inside intersection points on the line. If no inside points remain on the line, then the edge is ignored. The sign of the area is determined from the slope of the line. If the horizontal component of the line equation is positive, the sign of the area is positive. A negative horizontal component will result in a negative area. Vertical edges will have zero area. The areas are summed to come up with the final Thiesson area of the data point.

The three dimensional algorithm used for finding the Thiesson volume for a polyhedron clipped by the window is an extension of the two dimensional version. It involves first defining candidate faces for the polyhedron. These are the coefficients of the planes that potentially the faces of the Thiesson polyhedron will be constructed. Assuming the data points have been tesselated to maintain the Delaunay criterion, candidate faces are defined as the perpendicular bisecting planes to the tetrahedron edges radiating from the data point. The planes comprising the bounding parallelepiped are also defined as candidate faces. Each candidate face is intersected to define the edges of the polyhedron faces.

At this point, the problem has essentially been reduced to multiple two dimensional cases, where each face is treated individually. The edges on each face are intersected to define the candidate vertices. Upon arriving at a list of candidate vertices, the vertices are culled by defining them as inside or outside of all candidate faces. The remaining vertices define the vertices of the Thiesson polyhedron. Having defined the topology of the polyhedron, equations 6 and 7 can be used to compute the final volume. The area of each face of the Voronoi volume can now be determined and combined to compute the final volume.

At this point, all the geometric information needed to compute the local coordinates for each data point is known. All that remains to define all the variables needed for the Natural Neighbor interpolant (equation 2)

5. Estimating Gradients

Natural neighbor interpolation also uses an estimate of the gradient at each data point to increase the continuity of the function. This involves preprocessing each data point and computing and storing values for dz/dx and dz/dy or df/dx, df/dy, df/dz, for 3D. Computing the gradient estimate, like the interpolation itself, involves sampling neighboring data points to estimate the slope of the function at the location of a data point. The question, however arises as to which neighboring points should be used when computing the estimate and how they should be weighted.

A limited number of methods for computing the gradient were experimented with. It appears that many more methods could have been tested and their accuracy and robustness determined. A more complete study of gradients and how they affect interpolants for natural neighbor interpolation will be left to future research. The methods considered are discussed briefly.

5.1 Local Coordinates

[Sibson,81] suggests a method using local coordinates. This method for computing gradients seems to be most in harmony with the spirit of natural neighbor interpolation. With this method, gradients are determined only from the neighbors sharing common triangles. A local coordinate, l_m is defined for each neighboring vertex that defines its weight. First, Thiesson areas are computed for neighbor vertices to the data point in question. The point is then temporarily deleted from the TIN and new resulting Thiesson areas are computed. The local coordinate is defined exactly as is the local coordinate for the interpolation (see equation 1): Sibson defines the gradient vector bn as follows:

(9)

where

(10)

Although this method may seem to be the logical method to use, there appeared to some drawbacks. Computing the Thiesson area once the point has been deleted involves triangulating an arbitrary polygon to maintain Delaunay criterion. This polygon is not guaranteed to be convex. Although methods have been developed to solve just this problem, they are often time consuming. In addition, since the objective of this research was to develop a three dimensional algorithm, extending this case to three dimensions did not appear to be trivial as it would involve tessellating concave polyhedra into tetrahedron maintaining Delaunay criterion.

5.2 Neighboring Triangle Normals

One of the most straight-forward methods tried was to use the neighboring triangle normals to the data point. The normal at a data point is defined as the normalized sum of all of the triangle normals surrounding it. The gradient is then computed from the normal. In effect, this amounts to using the average gradient of the triangles surrounding the point. This method is very similar to how the Garoud shading model determines vertex normals.

Implementation of this method involved looping around the point where the gradient is to be estimated and computing the cross product of two edges of each adjacent triangle. Each component of the result is summed with the previous adjacent triangle normal. This is illustrated in equation 11.

(11)

where N is the estimated normal at the data point and n is the number of triangles adjacent to the data point. The gradient is then computed from the normal and stored with the vertex structure

Although this method was used for the 2D interpolation algorithm, this method also has some problems. An unnatural weighting to one side or other of the data point can occasionally occur if data points are not evenly distributed. This is most noticeable near the edges of the TIN.

Figure 9. Triangle group

Figure 9 shows a case where the estimated normal at the data point in the center would be unnaturally weighted to the right.

5.3 Least Squares Fit

This method involves estimating the gradient by using neighboring points to fit a plane through the data point in a least squares fashion. Neighboring data points are defined as those data points defining triangles whose circumcircle contain the point in question. The general formulation of the equations for the two dimensional case are shown in (12)

(12)

where Dx_i, Dy_i, Dz_i are the coordinate distances from the point in question to each neighboring data point. The gradients may be computed by solving a 2X2 system of linear equations.

To increase the accuracy of this method, a weight is also applied to each neighboring data point. The effect is to decrease the influence of a point as its distance increases. The weight, w of the neighboring point i is defined as:

(13)

where r is the greatest distance from the data point to any of its neighboring points and di is the distance to the neighboring point i. By including weights, the expression for the gradient can be refined to:

(14)

It appears that of the methods discussed, the final method is the most efficient and robust, particularly for cases illustrated in figure 9. The equivalent 3D equation for computing the gradient is shown here:

(15)

6. Extrapolation

One of the main advantages of natural neighbor interpolation is its ability to extend to extrapolation with little modification. The interpolated surface beyond the convex hull of the data points will maintain its first or second derivative continuity. In spite of these benefits there is a case worthy of note.

As was discussed previously, natural neighbor interpolation assigns weights based on the Thiesson area or volume of neighboring data points. The Thiesson polytopes for data points on the edge of the TIN must be clipped by the bounding window. Because of this, the area or volume of Thiesson polytopes at the boundary of the tesselation are greatly influenced by the size of the bounding window. In general, the larger the area of a Thiesson polygon the greater its relative weight on the interpolant. This means that the relative weight applied to an interpolant located beyond the convex hull from a data point on the edge of the TIN will change depending on the location of the boundary window.

Physically, this means that if a large slack area is provided between the outer data points and the boundary window, the weights defined at the data points on the TIN boundary will dominate as the size of the boundary window increases. Because the weight applied to an extrapolation point from neighboring data points is dependent on the size of the window, the appearance of the final interpolated surface will vary somewhat depending on the boundary window selected.

In spite of this characteristic inherent with extrapolation using natural neighbor interpolation, the results seem to conform with intuition. Although, in the ideal, we would like to define a unique interpolated surface, regardless of the size of the bounding window, the edge effects are an acceptable consequence of this method and perfectly acceptable for the prescribed application.

7. Results

Several test cases were run using the code generated from the 2D implementation of natural neighbor interpolation. The input for all test cases was a set of scattered x-y-z data points. One of the test cases is documented here. This test case is taken from actual seismic field data. It illustrates the advantages of natural neighbor interpolation when used for interpolating data arranged along lines. Figure 10a shows the triangulated data points used as input to the interpolation. Figures 10b and 10c compare two other common triangle-based interpolation methods. Both of these methods do not perform well when long-skinny triangles appear in the Delaunay triangulation. Additionally, the methods shown do not extend easily to extrapolation. Figure 10d shows contours from the same data points using natural neighbor interpolation. This figure illustrates the resulting smooth contours typical with natural neighbor interpolation. Both extrapolation and interpolation were used to compute figure 10d.

	Figure 10a. Triangulated data points used for example.
	Figure 10b. Contours based on linear triangle based interpolation.
	Figure 10c. Contours based on Clough-Tocher interpolation of the data points
	Figure 10d. Contours based on natural neighbor interpolation and extrapolation of the data points

There were several test cases generated to test the accuracy and efficiency of the three dimensional implementation of natural neighbor interpolation. For each of these cases, a set of three dimensional scattered data points along with function values was provided as input to the program. The objective of the program was to compute a specified number of interpolation or extrapolation points based on the input data.

Similar to the two dimensional implementation, the 3D interpolation program served as a preprocessor for a three dimensiona graphics package. Iso-parametric surfaces may were computed at specified function values or cross sections cut from the volume to expose interior data. Contours or fringes were also displayed on the surface of the cross sections. All of these methods were utilized in evaluating the interpolation method.

One of the cases tested samples a known algebraic function whose data is arranged along vertical traces. This arrangement of data is similar to what would be expected from borehole analysis. The control function is as follows:

(16)

The data is arranged in six lines of 10 data points each for a total of sixty data points. The Delaunay tesselation of the data points is shown in figure 11a.

	Figure 11a. Triangulated data points used for example.
	Figure 11b. Cut-away view of iso-surfaces of the control function.
	Figure 11c. Cut-away view of iso-surfaces based on natural neighbor interpolation of data points.
	Figure 11d. Cut-away view of iso-surfaces of the control function.
	Figure 11e. Cut-away view of iso-surfaces based on natural neighbor interpolation of data points.

8. Conclusion

It has become clear through this research that there currently exists a wide range of interpolation methods in both two and three dimensions. It is unlikely that any one of these methods will be able to accurately reproduce an underlying function in all circumstances. As was mentioned, there are several characteristics that are desirable in an interpolation scheme, some of which will change depending on the specific application.

Through this research we have been able to determine that natural neighbor interpolation has many properties which can make it a useful tool particularly in the field of site characterization. Some of these are listed following:

1. Generation of a C¹ or C² continuous function resulting in an esthetically pleasing surface or volume.
2. The local nature of the interpolation, allowing for relatively constant times to compute each interpolant, regardless of the number of data points in the data set.
3. The ability to extrapolate beyond the convex hull of the data points.
4. The ability to take into account both geometric and topological factors of the arrangement of the data points. This permits it to more accurately reproduce a function when unusual arrangement of data points is encountered.

It is clear that natural neighbor interpolation will not solve all problems, however this research has demonstrated some of the unique advantages of this interpolation technique. For at least the three dimensional case, the authors were unable to discover any published information. In order to fully implement the concepts first introduced by Sibson for the two dimensional case, new algorithms and techniques had to be developed to extend them to three dimensions.. Some important ground work has been laid and a useful scientific visualization tool developed.

The computer code developed in conjunction with this research, although adequate in many respects could be enhanced to speed up execution times. It appears that the majority of the time needed to compute the natural neighbor interpolant is taken in computing the volume of neighboring Thiesson polyhedra. For a typical data set, literally millions of manifold intersections must be computed. Methods should be determined where redundant calculations or simplifications can be made.

The properties of natural neighbor interpolation make it ideal for application in the field of geotechnical engineering. The computer code developed as part of this research could very easily take the place of many existing interpolation schemes. With enhancements added to speed up the code, the times needed to reproduce an underlying function from a data set could be dramatically reduced.

As a tool for geotechnical engineering, natural neighbor interpolation could also be customized to incorporate some of the specific requirements of this field. Since interpolation assumes that the underlying function is continuous, features such as faults and changes in stratigraphy are not taken into account. In addition, arbitrary shapes for interpolation regions may be necessary in order to take advantage of known boundary conditions such as bedrock or inclined or warped stratigraphy.

As was mentioned in the Introduction, three dimensional interpolation may be used to reproduce any continuous function in three dimensional space. Because of this the applications for three dimensional natural neighbor interpolation are numerous. The research presented here provides a potential contribution to many diverse fields.

References

Akima, Heroshi, A Method of Bivariate Interpolation and Smooth Surface Fitting for Irregularly Distributed Data Points, ACM Transactions on Mathematical Software, Vol. 4, No. 2, 1978, pp. 148-159.

Alfeld, Peter, A Trivariate Clough Toucher Scheme for Tetrahedral Data, Computer Aided Geometric Design, Vol. 1 1984, pp. 5-16

Alfeld, Peter, Scattered Data Interpolation in Three or More Variables, Mathematical Methods in Computer Aided Geometric Design, Tom Lynche and Larry L. Schumaker eds., Academic Press, Boston 1989, pp. 1-33

Baker, T. J., Automatic Mesh Generation For Complex Three-Dimensional Regions Using A Constrained Delaunay Triangulation, Engineering With Computers, Springer-Verlag New York, Vol. %, No. 3/4, Summer/Fall 1989, pp. 161-75

Bowyer, A., Computing Dirichlet Tessellations, The Computer Journal, Vol. 24, No. 2 1981, pp. 162-66

Davis, J. C., Statistics and Data Analysis in Geology, John Wiley & Sons, New York, 1986

Dahman, Wolfgang, Smooth Piecewise Quadric Surfaces, Mathematical Methods in Computer Aided Geometric Design, Tom Lynche and Larry L. Schumaker eds., Academic Press, Boston, 1989, pp. 181-193.

De Floriani, Leila, Bianca Falcidieno, and Caterina Pienovi, Delaunay-Based Representation of Surfaces Defined Over Arbitrary Shaped Domains, Computer Vision, Graphics, and Image Processing, Vol. 32, pp. 127-140.

Dyn, Nira, David Levin and Samuel Rippa, Data Dependant Triangulations for Piecewise Linear Interpolation, IMA Journal of Numerical Analysis, Vol. 10, 1990, pp. 137-154.

Franke, Richard, Scattered Data Interpolation: Tests of Some Methods, Mathematics of Computation, Vol. 38, No. 157, 1982, pp. 181-200

Franke, Richard and Greg Nielson, Smooth Interpolation of Large Sets of Scattered Data, International Journal for Numerical Methods in Engineering, Vol 15, 1980, pp 1691-1704

Goldman, Ronald N., Area of Planar Polygons and Volume of Polyhedra, Graphics Gems II, ed. James Arvo, Academic Press, New York, 1991.

Jones, Norman L., Applications of Computer-Aided Design Methods for Site Characterization in Civil Engineering, Thesis Geotechnical Engineering Center, Civil Engineering Department, University of Texas at Austin, May 1988

Jones, Norman L., Solid Modelling of Earth Masses for Applications in Geotechnical Engineering, Dissertation University of Texas at Austin, December 1990

Jones, T. A., D. E. Hamilton, and C. R. Johnson, Contouring Geologic Surfaces With the Computer, Van Nostrand Co.., New York, 1986

Lancaster, Peter and Kestutis Salkauskas, Curve and Surface Fitting, Academic Press, London, 1986.

Lawson, C. L., Software for C1 Surface Interpolation, Mathematical Software III, John R. Rice ed., Academic Press, New York 1977, pp. 161-94

Lee, D. T., and B. J. Schacter, Two Algorithms for Constructing A Delaunay Triangulation, International Journal of Computer and Information Sciences, Plenum Press, New York, London, Vol. 9, No. 3 June 1980, pp. 219-42

Petrie, G. and T. J. M. Kennie, Terrain Modeling in Surveying and Civil Engineering, Computer-Aided Design, Vol. 19, No. 4, 1987, pp. 171-187.

Rescorla, K. L., C1 Trivariate Polynomial Interpolation, Computer Aided Geometric Design, Vol. 4 1987, pp. 237-44

Sabin, M. A., Contouring - The State of the Art, Fundamental Algorithms for Computer Graphics, R. A. Earnshaw ed., Springer-Verlag, Berlin, Heidelberg, 1985, pp. 411-482.

Sibson, R., A Brief Description of Natural Neighbor Interpolation, Interpreting Multivariate Data, 1981, pp. 21-36

Watson, D. F., Computing the n-Dimensional Delaunay Tessellation with Application to Voronoi Polytopes, The Computer Journal, Vol. 8, No. 2 1981, pp. 167-72

Watson, D. F. and G. M. Phillip, Triangle Based Interpolation, Mathematical Geology, Vol. 16, No. 8 1984, pp. 779-95

Watson, D. F. and G. M. Phillip, Systematic Triangulations, Computer Vision, Graphics, and Image Processing, Vol. 26, No. 2, pp. 217-23

Watson, D. F. and G. M. Phillip, A Refinement of Inverse Distance Weighted Interpolation, Geo-Processing, Vol. 2, No. 4 1985, pp. 315-27

Watson, D. F. and G. M. Phillip, Neighbor Based Interoplation, Geobyte, Vol. 2, No. 2 1987, pp. 12-16