APPLICATIONS OF VORONOI AND DELAUNAY DIAGRAMS IN THE SOLUTION OF THE GEODETIC BOUNDARY VALUE PROBLEM

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.


INTRODUCTION
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., spacefrequency 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 andEscobar (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).

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).

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 m n > , 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 i π 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.  where i g Δ is the Helmert anomaly, the free-air anomaly corrected for the topographic and atmospheric attraction effects (e.g., MORITZ, 1984), and i a 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., ψ = 0), leading the indetermination of S(ψ).To avoid the indetermination, in this case, the Eq. ( 4) is replaced by (HEISKANEN & MORITZ, 1967), where Δ g 0 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 ) , ( λ φ P is split into three components (e.g., STRANG VAN HEES, 1986;DENKER and WENZEL, 1987;RAPP and BASÍC, 1992) ) , ( ) , ( ) , ( ) , (    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 ) , ( λ φ STOKES N .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.

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.5arcmin 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.

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 s i is the Euclidean distance from computation point P, with elevation H P , to generic integration cells i, with elevation H i and area a i .
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.
Figure 14 -Indirect effect for the Helmert's second condensation method, with the curve of -0.05 meter in highlight, as it was computed using Delaunay triangulation.

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, ∂Δg/∂H, 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 Δg P is referred to a fixed point P on which ∂Δg/∂r 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 R 2 dσ, expressed in terms of the angular distance ψ by 2 ψ 2Rsin l = . (10) 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 a i and mean gravity anomaly value Δ g i .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.
Figure 16 -Vertical gradient of Helmert gravity anomaly with the curve of zero mGal/km in highlight, as it was computed using Delaunay triangulation.

CONCLUSIONS
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.

Figure 1 -
Figure 1 -Voronoi (a) and Delaunay (b) diagram construction. Fig of long-wavelength components computed using geopotential model coefficients (e.g.,RAPP et al., 1991).The primary indirect effect term, 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.Both Delaunay and Voronoi schemes were used to compute 7 and 8 depict the results for the area of the Rio de Janeiro State.

Figure 10 -
Figure 10 -Geometry of a triangular prism with tilted top.

Figure 11 -
Figure 11 -Contribution per distance range in the terrain correction at the point of maximum value, up to 24 km from the point.

Figure 12 -
Figure 12 -Terrain correction on State of Rio de Janeiro area, with the curve of 2 mGal in highlight, as it was computed using Delaunay triangulation.

Figure 13 -
Figure 13 -Contribution per distance range in the indirect effect correction at the point of minumum value, up to 24 km from the point.

Figure 15 -
Figure 15 -Contribution per distance range in the vertical gradient of Helmert gravity anomaly at the point of minumum value, up to 24 km from the point.

Table 1 -
Differences (metres)between the classical numerical integration technique and Voronoi Delaunay.