Parallel Algorithm for Natural Neighbor Interpolation Alexander Tsidaev 1 Bulashevich Institute of Geophysics, Yekaterinburg, Russia 2 Ural Federal University, Yekaterinburg, Russia pdc@tsidaev.ru Abstract. This paper describes parallelization technique for Natural Neighbor interpolation algorithm. It is based on Green-Sibson Voronoi tesselation method and works without use of Delaunay triangulation. Obtained results show that algorithm is efficient enough to be used for sequential interpolation of multilayer data. Keywords: interpolation · natural neighbor · parallel algorithm · geo- physics · CUDA · OpenMP 1 Introduction Interpolation is very important part of many scientific processes, which rely on measured data. Obtaining of such a data is always complicated with different factors: usually it is lack of resources and physical impossibility to perform mea- surements in some point. Also, most of the interpretation methods work with regular grids while real data is usually measured as a set with irregular structure. This is why the interpolation is frequently needed to convert measured sets to some convenient initial models, which will be then interpreted. Particularly, author experiences the need of fast interpolation schemes during conversion of seismic profile data to 3D model on regular grids [1]. Example is presented on Fig. 1. Empirically it has been found, that the best results are obtained using Nat- ural Neighbor algorithm [2]. But on relatively large grids (1M points and more) this algorithm, which requires Voronoi diagram calculation for each point, works with insufficient efficacy. I.e. in 3D layered models the interpolation should be performed independently in a big set of layers (e.g. Fig. 1 consists of 801 layers: 80 km depth per 100 m). Even algorithms implemented in commercial software (such as Golden Software Surfer) took hours and even days to complete inter- polation for all layers. So the efficient parallel natural neighbor algorithm is required. Current paper presents the needed theory for the problem (Sections 2 and 3), algorithm itself (Section 4) and the result of the test in comparison with non-optimized implementation (Section 5). Parallel Algorithm for Natural Neighbor Interpolation 79 Fig. 1. Example of interpolation. Spare space between seismic profiles (a) is filled with interpolated values (b) 2 Voronoi Diagrams Let X be the 2-dimensional plane and {Pj } ∈ X be the set of the points (called seeds, sites or generators). Voronoi tesselation [3] of the X is a set of regions (polygons), each containing one seed and all points closer to that seed than to any other (Fig. 2): Rk = {x ∈ X | d(x, Pk ) ≤ d(x, Pj )} (1) Here Rk - element of Voronoi diagram; d(x, a) is a distance between points x and a and d(x, A) = inf{d(x, a) | a ∈ A}. Fig. 2. Voronoi diagram for 13 points (in white). Each element is presented in a different color It is easy to notice that all segments that form the boundaries between regions are the perpendiculars bisectors to the segments, which connect points of the initial set. 80 Alexander Tsidaev 3 Natural Neighbor Algorithm Natural Neighbor algorithm had been developed in 1981 by Robin Sibson [2]. Input data for the algorithm is a set of points {(xi , yi )}N i=0 given with some function f values, specified in these points: {f (xi , yi )}N i=0 . Algorithm idea is to calculate Voronoi diagram for all initial points (xi , yi ), and then to add each interpolated point (x, y) into the tessellation with sequential diagram recalcula- tion. The value G(x, y), which is attributed to the interpolated point, depends on how much of the area of initial diagram elements was “stolen” by the region of new inserted point (Fig. 3). N X G(x, y) = wi f (xi , yi ) (2) i=1 here f (xi , yi ) is the measured value in point (xi , yi ), wi = Q Rk is ratio of “stolen” k area (see Fig. 3). Rk is the area of the initial Voronoi diagram element for point Pk ; Qk is the intersection area of Rk and newly constructed element for the point (x, y). Fig. 3. Natural neighbor interpolation. The coloured circles, which represent the inter- polating weights, wi , are generated using the ratio of the shaded area to the cell area of the surrounding points. The shaded area is new Voronoi element for the point to be interpolated (based on image by Markluffel, distributed under CC BY-SA 3.0) So the natural neighbor algorithm is basically the algorithm to insert an additional point into existing Voronoi diagram. And it is well known that Voronoi diagram is a dual graph of Delaunay triangulation. This is why most of methods of the iterative Voronoi diagram construction are based on triangulation (e.g. [4], [5]). This approach has many advantages in case when many incremental steps should be performed sequentially. But in our case we always need to insert only one new point and to calculate resulting weights wi then. After this the obtained modified Voronoi diagram is no longer needed and we construct next Parallel Algorithm for Natural Neighbor Interpolation 81 one for the next point. The weight wi calculation depends on geometrical features of Voronoi diagram and could not be performed on its dual graph. This means that natural neighbor method requires reconstruction of Voronoi from Delaunay on each iteration step if implemented using the triangulation. In this research we use an algorithm for incremental Voronoi diagram construction by P.J. Green and R. Sibson [6] with natural neighbor weights calculation added. It does not require costly conversion of Voronoi graph to Delaunay graph and back on any step. 4 Iterative Algorithm for Voronoi Diagram Construction As described in previous section, the natural neighbor algorithm requires inser- tion of new point and recalculation of Voronoi diagram. To speed up the process and make easier understanding of how the existing polygons shapes were mod- ified, it is better to implement incremental point insertion algorithm instead of whole diagram recalculation. Let Pi be the initial points set with values f (Pi ). Let Ri be the Voronoi polygon for site Pi . We need to add new site P into the Voronoi diagram for Pi . Let d(Pa , Pb ) be the distance between points Pa and Pb . Then the algorithm can be written as follows. Algorithm 1 Incremental Voronoi diagram construction step of Natural Neighbor interpolation method 1: Select . Get the closest to P point Pn : d(P, Pn ) = min{d(P, Pi )} ∀i 2: Calculate perpendicular bisector to segment P Pn 3: Find two edges of Rn , which cross the bisector: E1 and E2 . Assume for definiteness that rotation direction of minimum angle between E1 and E2 with respect to point P is counter-clockwise. 4: Construct first edge of new Voronoi diagram element R for point P by connecting together points of intersection with bisector 5: Calculate coefficient wn as ratio of areas as described in previous section . Now constructing remaining edges of region R in counter-clockwise direction 6: finalEdge ← E1 7: while E2 6= finalEdge do . While rotation is not finished and we have not returned to the Rn 8: Find region Rk , which shares edge E2 with Rn 9: (E1 , E2 ) ← findIntersectedEdges(P, Rk ) . Same as steps 2-3 of the algorithm 10: Construct next edge of new region . See step 4 11: Calculate wk . See step 5 82 Alexander Tsidaev 12: end while PN 13: return i=1 wi f (xi , yi ) . All region edges were constructed and result value can be calculated 5 Parallel Implementation and Test Results Parallelization idea is quite simple and follows from the non-obstructive nature of Algorithm 1. Since the value in each interpolation point is calculated indepen- dently and none of the initial data is modified during process, then any number of point values can be calculated at the same time. Algorithm 1 was implemented in C++ for two parallel platforms: OpenMP and CUDA. Initial interpolation point set consisted of 644 points. Interpolation grid was a square of 1024x1024 points. All tests were performed on “Uran” super- computer (Krasovskii Institute of Mathematics and Mechanics, Yekaterinburg, Russia). Table 1. Interpolation time for a single layer (of 801) Hardware Calculation time, s Speed up factor 1 CPU core (Intel Xeon E5520) 187 1x 8 CPU cores (Intel Xeon E5520 OpenMP) 119 1.57x 1 GPU (Tesla M2050) 29.916 6.25x 2 GPUs (Tesla M2050) 17.096 10.94x 4 GPUs (Tesla M2050) 14.643 12.77x As it is seen from the Table 1, OpenMP parallelization showed unsatisfactory results. This may be connected with the fact that no attempts to perform manual optimizations were implemented and only the main loop was parallelized. Also it is important to note that other Intel-based optimization techniques, such as vectorization, are not applicable for current algorithm at all (because of branch- ing and different length of loops for different regions) - and so there is no chance to get additional speed up by introducing SIMD instructions. Acceleration obtained with GPU usage is quite acceptable, speed up factor is around one order of magnitude for current grid size and 4 GPUs. Our previous research showed that calculation time can have almost linear dependency versus number of GPU [7]. But in current case increase of involved GPUs does not en- tail corresponding linear increase of acceleration factor. This may be caused by complex program structure, which cannot be efficiently parallelized, and over- head connected with need to copy all the data to all GPUs. Program complexity caused the allocation of maximum allowed number of GPU registers, and this significantly reduced number of blocks that can be executed in parallel. Parallel Algorithm for Natural Neighbor Interpolation 83 6 Conclusion Parallel algorithm for natural neighbor interpolation was presented in this paper. Despite the obtained acceleration is enough for the current task of geophysical modeling, it seems perspective to implement algorithm for different paralleliza- tion platforms, such as MPI. Test example showed that even such complex al- gorithms, which cannot be vectorized well, can still be efficiently parallelized. And even MPI seems to be the most appropriate solution for parallel execution of inhomogeneous blocks, GPU computations could be used as a much cheaper alternative. References 1. Ladovsky I.V., Martyshko P.S., Druzhinin V.S., Byzov D.D., Tsidaev A.G., Kol- mogorova V.V.: Methods and results of crust and upper mantle volume density modeling for deep structure of the Middle Urals region. Ural Geophysical Messen- ger. 2(22). 31–45 (2013) 2. Sibson, R.: A brief description of natural neighbor interpolation (Chapter 2). In V. Barnett. Interpreting Multivariate Data. Chichester: John Wiley. 21–36 (1981) 3. Voronoi, G.: Nouvelles applications des paramtres continus la thorie des formes quadratiques. Journal fr die Reine und Angewandte Mathematik. 133 (133): 97–178 (1908) 4. Fan, Q., Efrat, A., Koltun, V., Krishnan, S., Venkatasubramanian, S.: Hardware- assisted natural neighbor interpolation. In C. Demetrescu, R. Sedgewick, R. Tamas- sia (Eds.), Proceedings of the Seventh Workshop on Algorithm Engineering and Ex- periments and the Second Workshop on Analytic Algorithms and Combinatorics. 111–120 (2005) 5. Guibas, L.J., Knuth, D.E., Sharir, M.: Randomized Incremental Construction of Delaunay and Voronoi Diagrams. Algorithmica 7. 381–413 (1992) 6. Green, P. J., Sibson, R.: Computing Dirichlet Tessellations in the Plane. The Com- puter Journal 21 (2): 168–173 (1978) 7. Tsidaev, A.: CUDA Parallel Algorithms for Forward and Inverse Structural Gravity Problems. Proceedings of the 1st Ural Workshop on Parallel, Distributed, and Cloud Computing for Young Scientists. 50–56 (2015)