On-line version ISSN 1982-2170
Bol. Ciênc. Geod. vol.18 no.3 Curitiba July/Sept. 2012
Applications of Voronoi and Delaunay Diagrams in the solution of the geodetic boundary value problem
Aplicações de Diagramas de Voronoi e Delaunay na solução do problema geodésico do valor de contorno
C. A. B. QuinteroI; I. P. EscobarII; C. F. Ponte-NetoI
IObservatório Nacional, Ministério da Ciência e Tecnologia, Rio de Janeiro, Brazil email@example.com; firstname.lastname@example.org
IIDepartamento de Engenharia Cartográfica, Universidade do Estado do Rio de Janeiro, Rio de Janeiro, Brazil email@example.com
Voronoi and Delaunay structures are presented as discretization tools to be used in numerical surface integration aiming the computation of geodetic problems solutions, when under the integral there is a non-analytical function (e. g., gravity anomaly and height). In the Voronoi approach, the target area is partitioned into polygons which contain the observed point and no interpolation is necessary, only the original data is used. In the Delaunay approach, the observed points are vertices of triangular cells and the value for a cell is interpolated for its barycenter. If the amount and distribution of the observed points are adequate, gridding operation is not required and the numerical surface integration is carried out by point-wise. Even when the amount and distribution of the observed points are not enough, the structures of Voronoi and Delaunay can combine grid with observed points in order to preserve the integrity of the original information. Both schemes are applied to the computation of the Stokes' integral, the terrain correction, the indirect effect and the gradient of the gravity anomaly, in the State of Rio de Janeiro, Brazil area.
Keywords: 2-D tessellation; Delaunay Triangulation; Voronoi Cells; Geodesy; Stokes' Integral
Este trabalho apresenta as estruturas de Voronoi e Delaunay como ferramentas de discretização a serem usadas na integração numérica de superfícies com o objetivo de resolver problemas computacionais geodésicos, quando no integrando a função não é analítica. No enfoque de Voronoi, a região de trabalho é particionada em polígonos, os quais contêm um ponto por polígono o que faz desnecessária a interpolação. No enfoque de Delaunay, os pontos observados são os vértices de um triangulo, e o valor da célula é o resultado da interpolação sobre o triangulo pelo seu baricentro. Se a quantidade de pontos observados e a sua distribuição são adequadas, a interpolação em grade não é necessária, e a integração é levada a cabo ponto a ponto. Mesmo quando a quantidade e distribuição dos pontos observados não são suficientes, as estruturas de Voronoi e Delaunay podem combinar grade com pontos observados a fim de preservar a integridade da informação original. Ambos os enfoques são aplicados ao cálculo da integral de Stokes, da correção de terreno, do efeito indireto, e do gradiente da anomalia da gravidade na região do estado de Rio de Janeiro, no Brasil.
Palavras-chave: Tesselação 2-D; Triangulação de Delaunay; Celulas de Voronoi; Geodésia; Integral de Stokes
In spite of the known techniques for geodetic parameters computing, either in the space-domain (SANTOS and ESCOBAR, 2000) or in frequency-domain (HAAGMANS et al., 1993; SCHWARZ et al., 1990), usually the target area is required to be partitioned into geographical grid elements. Geoid undulations, terrain corrections, indirect effects, for instance, are computed at these cells, based on gravity anomalies and heights, which are not evenly distributed. These are interpolated in order to produce a regular grid (HIRVONEN, 1956). Similarly, either integral methods (LEHMANN, 1997) or the combined ones (e.g., space-frequency domains), also require evenly distributed data (LI and SIDERIS, 1992). Gridding is usually (see Sideris 1995) required for fast Fourier transform (FFT) geoid determination techniques (e.g., STRANG VAN HEES, 1986; SIDERIS and TZIAVOS, 1988; HAAGMANS et al., 1993; FORSBERG and SIDERIS, 1993; KUROISHI, 2001), as well as for the fast Hartley transform (FHT) and fast T transform (FTT) techniques (AYHAN, 1997; LI and SIDERIS, 1992). Similarly, some traditional space-domain techniques, such as discrete summation (e.g., HIRVONEN, 1956; GEMAEL, 1999) and integral methods (e.g., JIANG and DUQUENNE, 1997; LEHMANN, 1997; NOVÁK et al., 2001), require the observed data to be gridded. Although this is not just a problem, nevertheless modified data are used instead of the original ones. Also, gridding usually expends excessive manual/computational effort, as well as it is liable to produce spurious data with the loss of genuine information (GOLDEN SOFTWARE, 1999). Besides, the interpolated value depends on the chosen gridding technique and on the grid 'nodes' separation, which are inherent to the spatial data distribution (e.g., GOLDEN SOFTWARE, 1999). In the worst case, it may produce spurious data that may lead to an inaccurate geoid. Also, different gridding methods provide different interpretations of the data because the grid node values are computed by different algorithms.
In this work, the discretization by means of Voronoi and Delaunay structures (AURENHAMMER, 1991; WATSON, 1981) is used to the computation of the Stokes' integral, terrain correction, indirect effect and gradient of the Helmert anomaly. Both approaches are reported in Santos and Escobar (2006, 2010). Genuine data are used and preserved if they have such a spatial distribution that does not require filling blank areas with interpolated data. In spite of a natural "smoothing" due to the data distribution, Voronoi and Delaunay schemes avoid the "synthetic" smoothing due to an interpolation step. Alternatively, if the original data are not sufficient, they can still be preserved in combination with a grid used to fill the blank areas. In Voronoi scheme, the target area is subdivided into a unique set of convex and adjacent polygonal cells, in which each one holds an original data point. In Delaunay scheme, the area is tessellated into contiguous triangular cells (triangulated irregular network - TIN). Mean values are locally interpolated for the barycenter (centroid) from respective triangles' vertices and remains enclosed to each cell (GOLDEN SOFTWARE, 1999).
2. THE TESSELLATION METHOD
Several authors have investigated Voronoi and Delaunay structures, their properties and applications in many different scientific fields. Those structures are used in the analysis of elements that have to be partitioned into contiguous domains called "natural neighbors". This term was coined by mathematicians, when noticing the frequent instances of those relationships found in the nature. Many geometric natural shapes tend to be organized into Voronoi-like figures, such as the formed surface layer of mud that has dried and contracted, or the tessellated bony shell of some turtles and tortoises.
In gravimetry, Delaunay triangulation has been used to model the topography for terrain corrections computation, in which the relief is represented by triangular prisms (RUPERT, 1988). Lehmann (1997) used a triangulation structure to model equilateral spherical triangles for the evaluation of geodetic surface integrals.
A Voronoi diagram is also referred to as the Dirichlet tessellation, and might be viewed as a geometric complement (a duality) to Delaunay triangulation (AURENHAMMER, 1991; TSAI, 1993). The polygon vertices are associated with Delaunay triangles by the same construction rule - the circumcircle criterion or the Delaunay criterion (TSAI, 1993).
2.1. Construction of Delaunay Triangulation and Voronoi Diagram
A Delaunay triangulation (also called a Delaunay simplicial complex) is a partition of an m-dimensional space, S, into adjacent triangular elements (Figure 1b). Given a set of n distinct points in S, so that n > m, every circumcircle associated with the triangles must contain no other points inside it. Whereas k points belonging to the perimeter of the area (peripheral points), the number of triangles generated is equal to 2(n - 1)- k. The interesting property of this structure (the approximated equiangular form) indicates that minimum angles are maximized and maximum angles are not minimized, which is an advantage over any triangulation of the same set of points. By Delaunay criteria, any triangulation with no obtuse angles has to be a Delaunay triangulation. According to Aurenhammer (1991), a triangulation without extreme angles (or "compact") is desirable, especially in methods of data interpolation. Additionally, in graphical computation the equiangular property is a need that provides the best visualization for displaying figures.
The vertices in Delaunay triangulation are all and only the n points of S, and its circumcenters are the actual Voronoi polygon vertices, or Voronoi polytopes (Figure 1a). The remaining circumcenters not satisfying the 'empty circumcircle criterion' are discarded from the construction.
The automatic contouring of the points is according to the triangulation algorithm by D.F. Watson (RUPERT, 1988), which was modified to include the Voronoi polygons' computation, in which the topological data structures set up the relations between data points, edges and Delaunay triangles. The algorithm implicitly ensures a closed bounding area perimeter (the convex hull), but it does not preserve its outer limits because this information is not required for the triangulation.
A Voronoi diagram is a partition of space S into n polygons, so that each one contains just one point i of the set. As a rule, given a polygon, a subset X of its vertices is closer to the inside point than any other point in S. Formally, we have
where d is the point-to-vertex Euclidean distance function.
The Voronoi structure is unique in a sense that each polygon edge is exactly halfway between each pair of sites in S. For this reason, the i - 1 half planes give rise to convex polygons, each of which have at most i - 1 edges (AURENHAMMER, 1991). Also, any point on an edge is equidistant to two sites, and any vertex is equidistant to each three sites, thus forming a polygonal partition of the region. By construction, no polygon can be empty, and as a consequence, space is partitioned into exactly i polygons.
The topological data structures for Voronoi diagram construction are almost the same as in Delaunay triangulation, but in Voronoi diagram the sequence of vertices and polygons edges is necessary to ensure the same area of computation as in the Delaunay triangulation. Additionally, a check for both Voronoi and Delaunay constructions is performed to verify if the sum of plane areas of the figures is the same.
3. APPLICATION OF VORONOI AND DELAUNAY TECHNIQUES
A comprehensive approach of Stokes' and Molodensky's formulations as well as the main Boundary Value Problem of Geodesy (BVPG) alternatives can be found in Guimarães and Blitzkow (2011). Voronoi and Delaunay diagrams were applied to compute the Stokes' integral for the local gravimetric geoid determination in the Rio de Janeiro State (and nearby regions), Brazil. The relief is very rugged, and varies between 0 - 2,600 m (Figure 2). In the same area Delaunay triangulation was applied to the computation of terrain correction, indirect effect and the gradient of Helmert gravity anomaly.
The dataset includes 1940 terrestrial gravity stations filled out with 491 Geosat 5-arcmin resolution gravity anomalies (Figure 3). The data on land are along some roads and a kriging interpolation was used to fill in most of the blank areas between the roads to a 5-arcmin resolution grid, amounting to 3973 data points.
Figures 4 and 5 depict the goal area tessellated according to the Delaunay and Voronoi schemes, respectively. The process removes clustered data inside a circle of radius 2000 m, in order to avoid rather irregular cells. These clusters accounted to 744 points for both the schemes, remaining 3229 data points, and the same amount of Voronoi cells were produced. Delaunay tessellation gave rise to 6278 triangular cells, whose vertices are the data points.
3.1 The Stokes' Integral Computation
In general, the geodesists deal with two distinct surfaces to represent the figure of the Earth. The determination of the distance between them is the main goal of the geodetic sciences. One is a theoretical surface - the reference ellipsoid - adopted as the planimetric referential for geodetic and geophysical applications. The other one is the geoid, which is the most important equipotential surface of the Earth's gravity field, the closest of the Earth's physical surface. It is a real surface, might be materialized, and may be approximately viewed as the mean sea surface, supposedly extended through the continents. The geoid is used as the altimetric referential for engineering applications.
Stokes' established the theoretical basis for the gravimetric determination of the geoid, considering the variation of gravity at different points on the Earth's surface. Stokes' formula solves the problem assuming a global and continuous data distribution over the Earth, and is given by (STOKES, 1849)
where is the gravity anomaly, R is the radius of a spherical Earth, and is normal gravity value on the ellipsoid surface. The spherical distance and azimuth are polar co-ordinates of a specific point referring to the gravity station, and is the spherical Stokes' function. Figure 6 outlines the geometrical relationship between polar co-ordinates and an elemental area on the spherical surface.
From Figure 1, we get the relationship
Substitution of Eq. (2) in Eq. (1) yields the geoidal undulation N as a function of the elemental area .
In practice, Stokes' integral is replaced by a discrete summation, due to the impossibility of an ideal data distribution. A method proposed by Hirvonen (1956) solves the discretization problem in order to determine the geoid undulations. It subdivides the studied area into a regular geographical grid, and each grid cell contains an interpolated mean gravity anomaly that represents that cell for the discrete evaluation. The geoidal undulation might be written as
where is the Helmert anomaly, the free-air anomaly corrected for the topographic and atmospheric attraction effects (e.g., MORITZ, 1984), and is an individual area corresponding to the i-th cell.
In Scheme Voronoi the actual values observed in the data points are assigned to polygons, and the calculation point coincides with a point of integration (i.e., y = 0), leading the indetermination of S(y). To avoid the indetermination, in this case, the Eq. (4) is replaced by (HEISKANEN & MORITZ, 1967),
where is the gravity anomaly at the point, and ψ0 is the mean spherical distance between the point and the respective edges of the Voronoi polygon.
In Delaunay scheme, as the computation is performed for the triangles' vertices (data points), and the value for the cell is computed for its barycenter, Stokes' function singularity has gone (SANTOS and ESCOBAR, 2010).
The geoidal height computation at point is split into three components (e.g., STRANG VAN HEES, 1986; DENKER and WENZEL, 1987; RAPP and BASÍC, 1992)
The represents the contribution of long-wavelength components computed using geopotential model coefficients (e.g., RAPP et al., 1991). The primary indirect effect term, , is computed by means of Helmert's second condensation method (Lambert 1930) from the elevation data file (Fig. 2). Finally, the component is computed from Eq. (4) with the residual free-air anomalies (free-air anomalies minus anomalies derived from model), corrected for the topographic relief and atmospheric attraction effects.
Involving almost the double of discretization cells, the Delaunay scheme provides a more smoothed aspect in component than the Voronoi scheme, what leads to a residual difference (Voronoi minus Delaunay) as is indicated in figure 9. The discrepancies range from -35 cm to 14 cm, with mean value of -2 cm and standard deviation of 4 cm.
A comparison between the application of Delaunay and Voronoi schemes and the classical technic in geoidal heights computation was done. As the first two terms in Eq. (6) are independent of the used method to compute Stokes' integral, because they do not depend on the way as the discretization is carried out, to examine the differences it is sufficient to consider only the term .
Table 1 presents the statistics of those differences for the Rio de Janeiro dataset. The RMS differences indicate an almost perfect agreement (99%) between the Voronoi and Delaunay schemes. Only one point in the Voronoi scheme is outside the 99% confidence interval for the RMS difference, which the maximum difference is less than 10 cm. With the exception of the area of Corcovado peak, results are statistically the same for the two methods.
3.2 The Terrain Correction Computation
The terrain correction takes into account the attraction effect of the irregularities of the topography in the vicinity of the gravity station and is always added to the observed value. Since the terrain correction can take values larger than other corrections to gravity (Earth's tide, free-air, Bouguer) it is very important, mainly in regions of rugged topography. The gravitational attraction of a vertical triangular prism is the mathematical model used here for the terrain correction computation (WOODWARD, 1975, RUPERT, 1988). The terrain is geometrically represented by vertical prisms of triangular base (B1, B2, B3) in horizontal plane at the altitude of the gravity point (P) and tilted top (A1, A2, A3) modeling the topography (Figure 10). Vertices of the triangular base are determined according to the Delaunay structures and the density is assumed to be constant. Within a given radial distance from each gravity point, the calculation of the vertical component of attraction of the prisms is computed.
Terrain corrections were computed for State of Rio de Janeiro area in a 1.5-arcmin grid, using Delaunay triangulation (Figure 4) and the dataset presented in Figure 3. A DTE dataset available online by the National Institute for Space Research, INPE, Brazil, was used to fill in the blank areas. This topographical database was unified and structured for the whole Brazil area - the TOPODATA Project (VALERIANO and ALBUQUERQUE, 2008). It was produced by a refined combination of local terrain elevation data and topographical information derived from Shuttle Radar Topography Mission - SRTM data (USGS, 2007), for a 1-arcsec horizontal resolution grid. A subset DTE grid of a 3-arcsec horizontal resolution was used to the radial distance of 3 km and a 15-arcsec resolution subset grid for radial distances between 3 km and 24 km. Distances and azimuths were calculated using Vincenty formulas for solving the inverse problem of Geodesy (VINCENTY, 1975). The values range between 0 and 37.22 mGal, with mean value of 1.66 mGal and standard deviation of 2.60 mGal. The Figure 11 presents a graphic with the contribution, in mGal, per distance range, in km, up to 24 km from the point of maximum value of terrain correction (37.22 mGal), it is possible observe that, for distances larger than 24 km, the contribution to terrain correction tends to zero. The Figure 12 shows the map of the terrain correction in Rio de Janeiro State area, computed using Delaunay triangulation (Figure 4) and the dataset presented in Figure 3.
3.3 The Indirect Effect Correction
As gravity reductions deal with topographic mass displacement, the resulting indirect effect on the geoidal height must be computed (HEISKANEN & MORITZ, 1967). The Helmert's second method of condensation involves generally small indirect effects, given by (SIDERIS and SHE, 1995):
where si is the Euclidean distance from computation point P, with elevation HP, to generic integration cells i, with elevation Hi and area ai.
Using the same process applied to the terrain correction, with the same Delaunay scheme, dataset and scanned distance, the indirect effect for Helmert's second method of condensation was computed. The values range between -0.424 m and 0 m. The mean value and standard deviation are, respectively, of -0.025 m and 0.032 m.
The Figure 13 presents a graphic with the contribution, in mm, per distance range, in km, up to 24 km from the point of minimum value of indirect effect correction (-0.424 m), it is possible observe that, for distances larger than 21 km, the contribution tends to zero. A map for the indirect effect is shown in the Figure 14.
3.4 The Gradient Of The Helmert Gravity Anomaly
For the solution of some geodetic problems, involving downward or upward continuation of the gravity, it is useful to know the gradient of the gravity anomaly, , which can be derived, assuming H and the radius vector, r, in the same direction (HEISKANEN and MORITZ, 1967, page 115, Eq. 2-217):
where is referred to a fixed point P on which is to be computed; R is the mean radius of the earth; l is the spatial distance between the fixed point P and the variable surface element R2 , expressed in terms of the angular distance by
The Eq. (24) can be written as
without error larger than 0.0006 %.
Since the gravity anomaly is not known as a continuous function, a numerical integration, based on Eq. (10), may be used for computing its vertical gradient, by discretizing the earth surface in cells, that is:
where n is the number of cells over the integrating area around the point P, i is the order of a generic cell with area ai and mean gravity anomaly value . As the summation decreases very rapidly with increasing distance, it is sufficient to extend the numerical integration over the immediate neighborhood of the point P.
Vertical gradients of the Helmert gravity anomaly were computed for State of Rio de Janeiro area in a 1.5-arcmin grid, using Delaunay triangulation (Figure 4) and the dataset presented in Figure 3. A grid of 1.5-arcmin horizontal resolution was used to fill in the blank areas. The neighboring area was scanned up to a radial distance of 24 km. The values of the gradient vary from -68.6 mGal/km to 27.2 mGal/km.
The Figure 15 presents a graphic with the contribution, in mGal/km, per distance range, in km, up to 24 km from the point of minimum value of vertical gradient (-68.6 mGal/km), it is possible observe that, for distances larger than 22 km, the contribution tends to zero. A map for the vertical gradient of the Helmert gravity anomaly is shown in the Figure 16.
Voronoi and Delaunay structures have been applied as alternative discretization tools to compute numerical surface integration in geodetic problems solutions, when under the integral there is a non-analytical function. Both schemes were used to computing the Stokes' integral, while terrain correction, indirect effect and vertical gradient of the Helmert gravity anomaly were computed using Delaunay triangulation.
The main advantage of those schemes is to perform the computation without an interpolation grid, when the amount and distribution of data points are enough. Even when this condition is not satisfied, it is possible to merge data points with a grid of interpolated data, used to fill in the blank areas. In both cases, the original data are preserved. However, when merged data are used, it is important to check the consistency between the interpolated data grid and the original data. Any bias or conflict should be eliminated a priori, to avoid artificial effects on results.
Both structures, of simple and efficient geometrical constructions, are useful for the tessellation of a site in order to evaluate the geoidal undulations by means of the Stokes' technique. The Voronoi approach uses less discretization cells than the Delaunay triangulation, nevertheless, both schemes leads to the same results, which are somewhat more efficient than the classical method.
Although a test with Voronoi scheme could have been done to the computation of terrain corrections, indirect effects and vertical gradients of the Helmert gravity anomaly, it was the Delaunay triangulation used here, having in sight the best fit of the triangles to rugged surfaces. For the sake of comparison, a test with Voronoi scheme could have been done.
AURENHAMMER F. Voronoi diagrams: A survey of a fundamental geometric data structure. ACM Computing Surveys 23(3): 345-405,1991. [ Links ]
AYHAN M.E. Updating and computing the geoid using two dimensional fast Hartley transform and fast T transform. Journal of Geodesy 71: 362-369,1997. [ Links ]
DENKER H., WENZEL H.G. Local geoid determination and comparison with GPS results. Bulletin Géodésique 61: 349-366,1987. [ Links ]
FORSBERG R., SIDERIS M.G. Geoid computations by the multi-band spherical FFT approach. Manuscripta Geodaetica 18: 82-90,1993. [ Links ]
GEMAEL C. Introdução à Geodésia Física. 1999. UFPR, Curitiba, Brazil [ Links ]
GOLDEN SOFTWARE Surfer Version 7 User Guide. 1999. Golden Software Inc., Golden, Colorado [ Links ]
GUIMARÃES G. do N., BLITZKOW D. Problema de valor de contorno da Geodésia: uma abordagem conceitual. Boletim de Ciências Geodésicas vol.17 no.4, p. 607-624, Curitiba Oct./Dec., 2011. [ Links ]
HAAGMANS R.R., DE MIN E., VAN GELDEREN M. Fast evaluation of convolution integrals on the sphere using 1D FFT, and a comparison with existing methods for Stokes' integral. Manuscripta Geodaetica 18: 227-241,1993. [ Links ]
HEISKANEN W.H., MORITZ H. Physical Geodesy. 1967. WH Freeman & Co, San Francisco [ Links ]
HIRVONEN R.A. On the precision of the gravimetric determination of the geoid. Transactions - American Geophysical Union 3 (1): 1-8, 1956. [ Links ]
JIANG Z., DUQUENNE H. On fast integration in geoid determination. Journal of Geodesy 71: 59-69, 1997. [ Links ]
KUROISHI Y. An improved gravimetric geoid for Japan, JGEOID98, and relationships to marine gravity data. Journal of Geodesy 74 (11-12): 745-765, 2001. [ Links ]
LAMBERT W.D. The reduction of observed values of gravity to sea level. Bulletin Géodésique 26: 107-181, 1930. [ Links ]
LEHMANN R. Fast space-domain evaluation of geodetic surface integrals. Journal of Geodesy 71: 533-540, 1997. [ Links ]
LI J., SIDERIS M.G. The fast Hartley transform and its application in physical geodesy. Manuscripta Geodaetica 17: 381-387, 1992. [ Links ]
MORITZ H. Geodetic Reference System 1980, Bulletin Géodésique 58 (3): 388-398, 1984. [ Links ]
RAPP R.H., WANG Y.M., PAVLIS N.K. The Ohio State University 1991 geopotential and sea surface topography harmonic coefficient models. 1991. Rep 410, Department of Geodetic Science and Surveying, The Ohio State Univ, Columbus [ Links ]
RAPP R.H., BASÍC T. Ocean-wide gravity anomalies from GEOS-3, Seasat and Geosat altimeter data. Geophysical Research Letters 19(19): 1979-1982, 1992. [ Links ]
RUPERT J. A gravitational terrain correction program for IBM compatible personal computers. 1988.Vol. 2.21, Geological Survey of Canada, GSC, Open File 1834 [ Links ]
SANTOS N.P., ESCOBAR I.P. Gravimetric geoid determination in the municipality of Rio de Janeiro and nearby region. Brazilian Journal of Geophysics 18 (1): 49-62, 2000. [ Links ]
SANTOS N.P. and ESCOBAR I.P.,2010. Use of EGM08 model and Shuttle Radar Topography Mission data for geoid computation in the State of Rio de Janeiro, Brazil: a case of study with Voronoi/Delaunay discretizations. Studia Geophysica et Geodaetica, 54: 239-255 [ Links ]
SIDERIS M.G., TZIAVOS I.N. FFT-evaluation and applications of gravity-field convolution integrals with mean and point data. Bulletin Géodésique 62: 521-540, 1988. [ Links ]
SIDERIS M.G., SHE B.B. A new high-resolution geoid for Canada and part of the U. S. by the 1D-FFT method. Bulletin Géodésique 69: 92-108, 1995. [ Links ]
SIDERIS M.G. Fourier geoid determination with irregular data. Journal of Geodesy 70: 2-12, 1995. [ Links ]
STOKES G.G.) On the variation of gravity on the surface of the Earth. In: Mathematical and Physical Papers, Vol. II, New York, pp 131-171 (from the Transactions of the Cambridge Philosophical Society, Vol. VIII, pp 672-695) , 1849. [ Links ]
STRANG VAN HEES G.L. Precision of the geoid, computed from terrestrial gravity measurements. Manuscripta Geodaetica 11: 86-98, 1986. [ Links ]
TSAI V.J.D. (1993) Fast topological construction of Delaunay triangulations and Voronoi diagrams. Computers & Geosciences 19(10): 1463-1474 [ Links ]
VINCENTY T. Direct and inverse solutions of geodesics on the ellipsoid with application of nested equations. Survey Review, XXII(176),88 - 93, 1975. [ Links ]
WOODWARD D.J. The gravitational attraction of vertical triangular prisms, Geophysical Prospecting, 23(3): 526-532, 1975, 1975. [ Links ]
Recebido em abril de 2012
Aceito em julho de 2012