See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/228725404

Interpolation as a tool for the modelling of three-dimensional geoscientific

datasets

Article · January 2005

CITATIONS READS

5 373

2 authors, including:

Christopher M. Gold

130 PUBLICATIONS   2,452 CITATIONS   

SEE PROFILE

Some of the authors of this publication are also working on these related projects:

Industrial Research Chair in Geomatics for Forest Applications View project

All content following this page was uploaded by Christopher M. Gold on 20 May 2014.

The user has requested enhancement of the downloaded file.

 INTERPOLATION AS A TOOL FOR THE MODELLING OF THREE-DIMENSIONAL
GEOSCIENTIFIC DATASETS

H. Ledoux∗, C. Gold

GIS Research Centre, School of Computing, University of Glamorgan

Pontypridd CF37 1DL, Wales, UK.
hledoux@glam.ac.uk — christophergold@voronoi.com

KEY WORDS: Modelling, Three-dimensional, Algorithms, Visualisation.

ABSTRACT:

Although interpolation methods are commonly used in geographical information systems (GIS) for modelling two-dimensional fields,
their three-dimensional counterparts are yet to be found in commercial software, and this despite the fact that more and more 3D datasets
are available, especially in geosciences. In this paper, we discuss the generalisation to 3D of four weighted-average interpolation
methods that are commonly used for 2D fields. We evaluate their properties when they are used with geoscientific data, and we also
discuss implementation details. We restrict our analysis to methods that can be used in a context of interactive exploration of a dataset.

1 INTRODUCTION modelling (Anselin, 1999; Bailey and Gatrell, 1995; Gold, 1993).
This is an approach where the user does not only use standards
Interpolation methods are an important part in a geographical in- operations (e.g. GIS spatial analysis or statistical methods) on a
formation system (GIS) and have been used for years to model, dataset, but actually interacts with it by manipulating and trans-
among others, elevation data. They are crucial in the visualisation forming it and then looks at the consequences of his manipula-
process (e.g. generation of contours lines), for the conversion of tions. The emphasis is put on the interaction between the user and
data from one format to another (e.g. from scattered points to the dataset through modification and visualisation tools. A sim-
raster), to have a better understanding of a dataset or simply to ple example of this approach is a spreadsheet tool (e.g. Microsoft
identify ‘bad’ samples. The result of interpolation—usually a Excel) where numerical data can be explored with the help of
surface that represents the real terrain—must be as accurate as different statistical tests and charts (histograms, scatter points,
possible because it often forms the basis for spatial analysis, for pie charts, etc.). These charts are linked to the data, i.e. if some
example runoff modelling or visibility analysis. Although inter- numbers in the dataset are changed, the charts are automatically
polation helps in creating three-dimensional surfaces, it is intrin- updated. In the context of geoscientific modelling, one can think
sically a two-dimensional operation for only the x−y coordinates of an environment where the samples can be inserted, deleted or
of each sample are used and the elevation is considered as an at- even moved at will by the user, and he gains insight about the
tribute. spatial variability of the field studied by manipulating the dataset
and observing the results. A key factor in interactive analysis
With the new technologies available to collect information about is the speed at which each function is performed, to ensure that
the Earth, more and more three-dimensional data are collected. A the user gets an (almost) instantaneous result from a query or
typical dataset in geosciences has samples in three-dimensional an operation, so that he is not disturbed by long waits. We dis-
space (x − y − z coordinates) to which one or more attributes are cuss in Section 2 the properties that an interpolation method for
attached. The attribute collected is usually a scalar, and typical geoscientific data should have, and we present in Section 3 four
examples are the salinity or the temperature of the water, the per- methods. For each, we present and evaluate their properties in 3D
centage of gold in the rock, or the humidity of the air. Samples space, and we also discuss implementation details to ensure that
are actually collected to study the spatial variability of an attribute they are computationally efficient.
over a certain extent in space; if a is the attribute studied, a contin-
uous function a = f (x, y, z) is defined, this is called a field (see
2 WHAT IS A GOOD INTERPOLATION METHOD?
Goodchild (1992) for more details). Because of the way they are
collected, three-dimensional geoscientific datasets often have a
highly anisotropic distribution. Geologic and oceanographic data Given a set of samples to which an attribute a is attached, spatial
are, for example, respectively gathered from boreholes and water interpolation is the procedure used to estimate the value of the
columns: data are therefore usually abundant vertically but sparse attribute at an unsampled location x. To achieve that, it creates
horizontally. To model such datasets, three-dimensional interpo- a function f , called the interpolant, that tries to fit the samples
lation methods that consider the specificities of the data must be as well as possible. Interpolation is based on spatial autocorre-
used. While most of the interpolation methods used in GIS in- lation, that is the attribute of two points close together in space
tuitively extend to 3D, it is not obvious that they preserve their is more likely to be similar than that of two points far from each
properties or are appropriate for such datasets. other. Watson (1992), in his authoritative book, lists the essential
properties of an ‘ideal’ interpolation method for bivariate geosci-
In this paper, we discuss the generalisation to three dimensions entific datasets. These properties can realistically be present for
of some of the interpolation methods found in GIS or geosci- trivariate datasets, and are as follows:
entific modelling packages. We focus our attention on meth-
ods that can be applied in a context of dynamic or interactive
1. exact: the interpolant must ‘honour’ the data points, or ‘pass
∗ Corresponding author. through’ them.
 2. continuous: a single and unique value must be obtained at
each location. This is called a C 0 interpolant in mathemat-
ics.

3. smooth: it is desirable for some applications to have a func-

tion for which the first or second derivative is possible every-
where; such functions are respectively referred to as C 1 and
C 2 interpolants.

4. local: the interpolation function uses only some neighbour-

ing samples to estimate the value at a given location. This
ensures that a sample with a gross error will not propagate Figure 1: The VD for a set of points in the plane.
its error to the whole interpolant.

5. adaptability, i.e. the function should give realistic results

for anisotropic data distributions and/or for datasets where
the data density varies greatly from one location to another.

In the context of dynamic/interactive modelling, we add to this

list:

6. the method must be computationally efficient.

7. the method must require as little input as possible from the

user, i.e. it should be automatic and not rely on user-defined
parameters that require a priori knowledge of the dataset. Figure 2: Two 3D Voronoi cells adjacent to each other.

cell have the same value (see Figure 1 and Figure 2 for respec-
We describe in the next section four interpolation methods, com-
tively the 2D and 3D examples). The Voronoi cell of a point
monly found in GIS and geoscientific software, that respect some
p ∈ S, defined Vp , is the set of points x ∈ Rd that are closer to p
of the criteria of an ideal method for an interactive system. The
than to any other point in S. In 2D, Vp is a convex polygon, and in
major omission from our list is Kriging (Matheron, 1971; Oliver
3D it is a convex polyhedron. Observe that the size and the shape
and Webster, 1990), which, although possessing many of the ideal
of a Voronoi cell is determined by the distribution of the points:
properties and being popular in many application domains be-
when the data is dense the cells are smaller. The union of the
cause it produces an interpolant that minimises the error variance
Voronoi cells of all generating points p ∈ S form VD(S). The
at each location, cannot be realistically be used in an dynamic en-
VD actually creates a piecewise model, where the interpolation
vironment. Indeed, the Kriging estimation is based on a function
function inside each Voronoi cell is a constant function.
characterising the dependence between the attributes of any two
samples that are at a given distance from each other. This function Although the method possesses many of the properties listed in
is obtained by studying the variance of the attributes according to Section 2 (it is exact, local and can handle anisotropic data dis-
the distance and the direction, and fitting a simple function. This tributions), it cannot be used with geoscientific data if realistic
is done manually by the user, and is a somewhat time-consuming results are wanted because it fails lamentably properties 2 and 3.
process that is impossible to do every time the dataset is modified. The interpolation function is indeed discontinuous at the border
of cells. It can nevertheless be useful in some applications: it is
The four interpolation methods discussed are all weighted-avera- for example often used in remote sensing to avoid averaging or
ge methods, which are methods that use only some sample points, blurring the resulting image, and it is also useful for nominal or
to which weights (◦ importance) are assigned. The interpolation ordinal data, e.g. rock types.
function f of such methods have the following form:
P The implementation of the method sounds easy: simply find the
P
k
i=1 wi (x) ai closest data point and assign its value to the interpolation loca-
f (x) = k
(1)
i=1 wi (x)
tion. The difficulty lies in finding an efficient way to get the clos-
est data point. The simplest way consists of measuring the dis-
where wi (x) is the weight of each neighbour pi (with respect to tance for each of the n points in the dataset, but this is too slow
the interpolation location x) used in the interpolation process, and for large datasets. To speed up this brute-force algorithm, aux-
ai the attribute of pi . iliary data structures that will spatially index the points must be
used. They usually subdivide hierarchically the space into cells
(usually squares or rectangles) and build a tree that uses O(n)
3 FROM TWO TO THREE DIMENSIONS space; examples of such trees are the KD-tree, the R-tree and the
octree (see van Oosterom (1999) for a survey of two-dimensional
3.1 Nearest Neighbour structures: the three-dimensional methods are simple extensions).
Thus, to find the nearest neighbour, it suffices to navigate in the
Nearest neighbour is a simple interpolation method: the value of tree and simply test the points in the adjacent cells. Locating the
an attribute at location x is simply assumed to be equal to the cell containing a test point can for example be done efficiently (in
attribute of the nearest data point. Given a set S of data points O(log n) time) with a KD-tree.
in a d-dimensional space, if interpolation is performed with this
method at many locations close to each other, the result is the Another solution consists of building the VD for the set of points
Voronoi diagram (VD) of S, where all the points inside a Voronoi and identifying inside which cell the interpolation point lies. In
 p1 p1
p3
p2 A2 p3
p4 x
wi (x) = |xpi |−h A3
A1
where h defines the x
p1 x p3
influence of pi on x wi (x) = Ai
p2
and |ab| is the distance
p5 p2 p4
2D 3D
(a) Figure 4: Barycentric coordinates.

(b)
Figure 5: Delaunay triangulation of the same set of points as in
Figure 3: (a) Inverse distance to a power interpolation. (b) Prob-
Figure 1.
lems with anisotropic datasets.
The generalisation of this method to three dimensions is straight-
fact, the Delaunay triangulation (DT)—which is the dual of the
forward: a searching sphere with a given radius is used. The
VD—can also be constructed as the two structures are equiva-
same problems with the one-dimensionality of the method will
lent, i.e. that the knowledge of one structure implies the knowl-
be even worse because the search must be performed in one more
edge of the other. Many algorithms exist to construct a three-
dimension. The method has too many problems to be considered
dimensional VD/DT, see for example the incremental methods of
for geoscientific datasets: the interpolant is not guaranteed to be
Watson (1981) or Edelsbrunner and Shah (1996). More details
continuous, especially when the dataset has an anisotropic distri-
about the latter method are available in Section 3.4. For a set of
bution, and the criterion has to be selected carefully by the user.
n points in R3 these algorithms construct the VD/DT in O(n2 ),
which is worst-case optimal since the complexity of the VD/DT
The implementation problems are also similar to the ones en-
for some distribution of points is quadratic. Finally, identifying
countered with the previous method, and an auxiliary data struc-
the cell containing a query point is very efficient and is expected
ture must be used to avoid testing all the points in a dataset.
to take only O(n1/4 ) when the points are randomly distributed
(see Mücke et al. (1999) for the details).
3.3 Linear Interpolation in Triangles
3.2 Distance-based Methods
This method is popular for terrain modelling applications and is
Distance-based methods are probably the most known methods based on a triangulation of the data points. As is the case for the
and they are widely used in many fields. As shown in Figure 3a, VD, a triangulation is a piecewise subdivision of the space cov-
in two dimensions they often use a ‘searching circle’, whose ra- ered by the data points, and in the context of interpolation a linear
dius is user-defined, to select the data points pi involved in the function is assigned to each piece (each triangle). Interpolating
interpolation at location x. The weight assigned to each is typ- at location x means first finding inside which triangle x lies (see
ically based on the square of the distance from x to pi . Other Mücke et al. (1999) for an efficient method), and then the height
weights can also be used, see Watson (1992) for a discussion. is estimated by linear interpolation on the plane defined by the
The size of the radius of the searching circle influences greatly three vertices forming the triangle. This can be efficiently im-
the result of the interpolation: a very big radius means that the plemented by using barycentric coordinates, which are local co-
resulting surface will be smooth or ‘flattened’; on the other hand, ordinates defined within a triangle (see Figure 4 for details). To
a radius that is too small might have dramatic consequences if for obtain satisfactory results, this method is usually used in 2D with
example no data points are inside the circle. A good knowledge a Delaunay triangulation because, among all the possible triangu-
of the dataset is thus required to select this parameter. lations of a set of points in the plane, it maximizes the minimum
angle of each triangle. In other words it creates triangles that are
This method has many flaws when the data distribution is highly as equilateral as possible (see Figure 5). A Delaunay triangula-
anisotropic or varies greatly in one dataset because a fixed-radius tion is a triangulation for which the circle circumscribed to each
circle will not necessarily be appropriate everywhere in the data- triangle is empty, i.e. it does not contain in its interior any other
set. Figure 3b shows one example where one circle, when used point in the set. This ensures that the three vertices used in the
with a dataset extracted from contour lines, clearly gives erro- interpolation process will most likely be close to and around the
neous results at some locations. The major problem with the interpolation location.
method comes from the fact that the criterion, for both selecting
data points and assigning them a weight, is one-dimensional and The generalisation of this method to 3D is as follows: linear
therefore does not take into account the spatial distribution of the interpolation is performed within each tetrahedron of a 3D tri-
data points close to the interpolation location. A criterion based angulation. Note that the barycentric coordinates can be gener-
on areas and volumes for respectively bivariate and trivariate data alised to linearly interpolate inside a tetrahedron. Finding ‘good’
yields better results, as Sections 3.3 and 3.4 explain. tetrahedra is however more difficult than finding good triangles
 p1

p2
w1
w2
x wi (x) = stolen area
p3 w3 w6
w4
p6
(a) p4 w5

p5
c

d Figure 7: Natural neighbour coordinates in 2D for the interpola-

a
e tion point x. The shaded polygon is Vx .
b
3.4 Natural Neighbour Interpolation

It has been shown by different researchers (Gold, 1989; Sam-

bridge et al., 1995; Watson, 1992) that the natural neighbour in-
(b) terpolation method (Sibson, 1981) avoids most of the problems
the other methods have with anisotropic datasets. This is a method
Figure 6: (a) A four-sided convex polygon can be triangulated by
based on the Voronoi diagram for both selecting the data points
two different ways, but the Delaunay criterion guarantees that the
involved in the process, and assigning them a weight. It uses two
triangles created will be as equilateral as possible. (b) In three
VDs: one for the set of data points, and another one where a
dimensions, five vertices can be triangulated with either two or point x is inserted at the interpolation location. Let us describe
three tetrahedra. Although the triangulation at the bottom has the method for the two-dimensional case. The insertion of x lo-
two nicely shaped tetrahedra, they are not Delaunay (the point cally modifies a VD: the Voronoi cell Vx of x ‘steals’ some parts
e is inside the sphere defined by the points a, b, c and d). The of some Voronoi cells, as shown in Figure 7. This idea forms
triangulation at the top is Delaunay, but contains one very thin the basis of natural neighbour coordinates (Sibson, 1980, 1981),
tetrahedron spanned by the points a, b, d and e. which define quantitatively the amount Vx steals from each of its
natural neighbours. Let D be the VD(S), and D+ = D ∪ {x}.
The Voronoi cell of a point p in D is defined by Vp , and Vp+ is its
because the max-min property of Delaunay triangles does not cell in D+ . The natural neighbour coordinate of x with respect
generalise to 3D. A 3D DT can indeed contain some tetrahedra, to a point pi is
called slivers, whose four vertices are almost coplanar (see Fig-
ure 6); interpolation within such tetrahedra obviously does not Area(Vpi ∩ Vx+ )
yield good results. Note that slivers do not have two-dimensional wi (x) = (2)
Area(Vx+ )
counterparts.
where Area(Vpi ) represents the area of Vpi . For any x, the value
of wi (x) will always be between 0 and 1: 0 when pi is not a nat-
Despite the occasional presence of slivers, Delaunay tetrahedra ural neighbour of x, and 1 when x is exactly at the same loca-
still have a more desirable shape than arbitrary tetrahedra. Ra- tion as pi . A further important consideration is that the sum of
jan (1991) shows that the smallest sphere containing a Delaunay the areas stolen from each of the k natural neighbours is equal
tetrahedron is smaller than that of any other tetrahedron, i.e. that to Area(Vx+ ). Therefore, the higher the value of wi (x) is, the
Delaunay tetrahedra tend to have a round shape. Also, Cheng stronger is the ‘influence’ of pi on x. The natural neighbour co-
et al. (2000) show that a triangulation where all the tetrahedra are ordinates are influenced by both the distance from x to pi and
well-shaped can be obtained as follows: start with the DT, assign the spatial distribution of the pi around x. Note that the method
some importance to some data points forming slivers, and modify is valid in any dimensions, for a Voronoi diagram exists in any
locally some tetrahedra (with the help of flips, see next section) to dimensions. Hence, volumes are stolen in three dimensions, and
obtain a regular triangulation of the set of points (Edelsbrunner hyper-volumes in higher dimensions.
and Shah, 1996).
Surprisingly, although many authors present the properties and
advantages of the method, few discuss details concerning its im-
As mentioned in Section 3.1, there exist different algorithms to plementation, and this is probably why its use is not widespread
construct and navigate efficiently in a 3D DT, and the linear in- among the GIS community. The two-dimensional case is rela-
terpolation itself is done in constant time since only four data tively easy to implement as efficient algorithms for construct-
points are involved. When used on a triangulation having well- ing a VD/DT (Fortune, 1987; Guibas and Stolfi, 1985; Watson,
shaped tetrahedra, the linear interpolation method possesses all 1981) and deleting a point from it (Devillers, 2002; Mostafavi
the ingredients of an ideal method but one: the first (or second) et al., 2003) exist. In three dimensions, things get more compli-
derivative of the resulting function is not possible at the vertices, cated because the algorithms to modify a VD/DT and computing
edges and faces of the triangulation. Akima (1978) solved this volumes are more complex (e.g. deleting a single vertex from a
problem in two dimensions by using higher order functions in DT is for example rather tricky to implement, as Ledoux et al.
each piece of the triangulation, and the same could be done in (2005) show). To our knowledge, the only efficient algorithm
three dimensions. is due to Boissonnat and Cazals (2002) and consists of mimick-
 a 4 DISCUSSION

flip23 We have shown that most of the weighted-average methods do

c generalise to 3D, but sometimes new problems appear, for exam-
b d flip32
ple with triangulations in three-dimensions or with the computa-
tional efficiency. Also, the flaws present in some two-dimensional
e methods will often be amplified in three and higher dimensions.
Natural neighbour interpolation seems to be the method that is the
Figure 8: flip23 and flip32. most appropriate for modelling geoscientific dataset because it is
robust, entirely objective, and gives good results when the data
distribution is highly irregular. We have already implemented in
ing the insertion of the interpolation point x to avoid deleting it two and three dimensions the algorithm described in this paper
afterwards. The stolen volume is obtained by a somewhat com- and we hope our method will make it possible for the GIS com-
plicated method that ‘orders’ the vertices surrounding x and then munity to take advantage of natural neighbour interpolation for
decomposes the volume into tetrahedra. The implementation of modelling geoscientific data.
this method is rather intricate.
Acknowledgements

We have therefore proposed a simple alternative that is valid in Part of this work was done when the first author was at the Hong
any dimensions (Ledoux and Gold, 2004). Our algorithm is based Kong Polytechnic University, and we thank the support of the
on the well-known flipping algorithm (Edelsbrunner and Shah, Hong Kong’s Research Grants Council (Project PolyU 5068/00E).
1996) to insert a point in a d-dimensional Delaunay triangula-
tion, and exploits the duality between a DT and the VD. Only
minor modifications to Edelsbrunner and Shah’s algorithm are REFERENCES
needed to compute the natural neighbour coordinates. In brief,
a flip is a local topological operation that modifies the configu- Akima, H., 1978. A method of bivariate interpolation and smooth
ration of adjacent tetrahedra in a triangulation (flips are actually surface fitting for irregularly distributed data points. ACM Trans-
possible in any dimensions, but for the sake of simplicity we only actions on Mathematical Software, 4(2):pp. 148–159.
discuss the 3D case here). Consider a set S = {a, b, c, d, e} of
points in R3 , as shown in Figure 8. There are two ways to trian- Anselin, L., 1999. Interactive techniques and exploratory spatial
gulate S: either with two or three tetrahedra. In the first case, the data analysis. In P. Longley, M. Goodchild, D. Maguire, and
two tetrahedra share a face bcd, and in the latter case the three D. Rhind, eds., Geographical Information Systems, chapter 17,
tetrahedra all have a common edge ae. A flip23 transforms a pp. 253–266. John Wiley & Sons.
triangulation of two tetrahedra into another one containing three
tetrahedra; a flip32 is the inverse operation. Also, it has been Bailey, T. C. and Gatrell, A. C., 1995. Interactive spatial data
shown that a single point can be inserted in a DT by applying a analysis. Wiley, New York.
sequence of flips (Joe, 1991). Our idea was simply to ‘remem-
ber’ the sequence of flips used to insert a point, and reverse it to Boissonnat, J.-D. and Cazals, F., 2002. Smooth surface recon-
delete it. We also showed that only the volume of some faces of struction via natural neighbour interpolation of distance func-
a Voronoi cell need to be computed in this context. The resulting tions. Computational Geometry—Theory and Applications,
algorithm is efficient (its time complexity is the same as the one 22:pp. 185–203.
for inserting a single point in a VD/DT), and we believe it to be
considerably simpler to implement than other known methods, as Cheng, S.-W., Dey, T. K., Edelsbrunner, H., Facello, M. A.,
only an incremental insertion algorithm based on flips, with some and Teng, S.-H., 2000. Sliver exudation. Journal of the ACM,
minor modifications, is needed. 47(5):pp. 883–904.

Devillers, O., 2002. On deletion in Delaunay triangulations.

International Journal of Computational Geometry and Applica-
Although computationally more intensive than the other methods,
tions, 12(3):pp. 193–205.
natural neighbour interpolation possesses all the properties enu-
merated in Section 2, except that the first derivative is undefined
Edelsbrunner, H. and Shah, N. R., 1996. Incremental topological
at the data points. To obtain a differentiable function everywhere,
flipping works for regular triangulations. Algorithmica, 15:pp.
Sibson uses the weights defined in Equation 2 in a quadratic equa-
223–241.
tion where the gradient at x is considered. To our knowledge, this
method has not been used with success with real data and there- Fortune, S., 1987. A sweepline algorithm for Voronoi diagrams.
fore we do not use it. Other ways to remove the discontinuities Algorithmica, 2:pp. 153–174.
at the data points have been proposed: Watson (1992) explains
different methods to estimate the gradient at x and how to incor- Gold, C. M., 1989. Surface interpolation, spatial adjacency and
porate it in Equation 1; and Gold (1989) proposes to modify the GIS. In J. Raper, ed., Three Dimensional Applications in Ge-
weight of each pi with a simple hermitian polynomial so that, as ographic Information Systems, chapter 3, pp. 21–35. Taylor &
x approaches pi , the derivative of f (x) approaches 0. Francis.

Gold, C. M., 1993. Forestry spatial decision support system clas-

Modifying Equation 1 to obtain a continuous function can yield sification, and the ‘flight simulator’ approach. In Proceedings
very good results in some cases, but with some datasets the re- GIS’93: Eyes on the Future, pp. 797–802. Vancouver, Canada.
sulting surface can contain unwanted effects. Different datasets
require different methods and parameters, and, for this reason, Goodchild, M. F., 1992. Geographical data modeling. Computers
modifications should be applied with great care. & Geosciences, 18(4):pp. 410–408.
 Guibas, L. J. and Stolfi, J., 1985. Primitives for the manipu- Rajan, V. T., 1991. Optimality of the Delaunay triangulation in
lation of general subdivisions and the computation of Voronoi Rd . In Proceedings 7th Annual Symposium on Computational
diagrams. ACM Transactions on Graphics, 4:pp. 74–123. Geometry, pp. 357–363. ACM Press, North Conway, New Hamp-
shire, USA.
Joe, B., 1991. Construction of three-dimensional Delaunay trian-
gulations using local transformations. Computer Aided Geomet- Sambridge, M., Braun, J., and McQueen, H., 1995. Geophysical
ric Design, 8:pp. 123–142. parameterization and interpolation of irregular data using natural
neighbours. Geophysical Journal International, 122:pp. 837–
Ledoux, H. and Gold, C. M., 2004. An efficient natural neighbour
857.
interpolation algorithm for geoscientific modelling. In P. Fisher,
ed., Developments in Spatial Data Handling—11th International Sibson, R., 1980. A vector identity for the Dirichlet tesselation. In
Symposium on Spatial Data Handling, pp. 97–108. Springer. Mathematical Proceedings of the Cambridge Philosophical Soci-
Ledoux, H., Gold, C. M., and Baciu, G., 2005. Flipping to ety, 87, pp. 151–155.
robustly delete a vertex in a Delaunay tetrahedralization. In Sibson, R., 1981. A brief description of natural neighbour inter-
Proceedings International Conference on Computational Science polation. In V. Barnett, ed., Interpreting Multivariate Data, pp.
and its Applications — ICCSA 2005, LNCS 3480, pp. 737–747. 21–36. Wiley, New York, USA.
Springer-Verlag, Singapore.
van Oosterom, P., 1999. Spatial access methods. In P. Longley,
Matheron, G., 1971. The theory of regionalised variables
M. Goodchild, D. Maguire, and D. Rhind, eds., Geographical
and its applications. Les Cahiers du Centre de Morphologie
Information Systems, volume 1, chapter 27, pp. 385–400. John
Mathématique de Fontainebleu, 5.
Wiley & Sons.
Mostafavi, M. A., Gold, C. M., and Dakowicz, M., 2003. Delete
and insert operations in Voronoi/Delaunay methods and applica- Watson, D. F., 1981. Computing the n-dimensional Delaunay tes-
tions. Computers & Geosciences, 29(4):pp. 523–530. sellation with application to Voronoi polytopes. Computer Jour-
nal, 24(2):pp. 167–172.
Mücke, E. P., Saias, I., and Zhu, B., 1999. Fast randomized point
location without preprocessing in two- and three-dimensional Watson, D. F., 1992. Contouring: A Guide to the Analysis and
Delaunay triangulations. Computational Geometry—Theory and Display of Spatial Data. Pergamon Press, Oxford, UK.
Applications, 12:pp. 63–83.
Oliver, M. and Webster, R., 1990. Kriging: a method of interpola-
tion for geographical information systems. International Journal
of Geographical Information Systems, 4:pp. 313–332.

View publication stats