Research paper
Approximation of Voronoï cell attributes using parallel solution

https://doi.org/10.1016/j.advengsoft.2019.03.012Get rights and content

Highlights

  • Rasterized approximation of hyper-dimensional Voronoi diagram.

  • Exctraction of selected scalar topological attributes of cells.

  • Efficient massively parallel solution providing significant speedup in comparison with full analytic solution.

  • Explicit estimation of rasterization error in general dimension.

Abstract

This paper concerns an algorithm for fast parallel approximation of selected attributes of hyper-dimensional Voronoï cells in a unit hypercube. The presented algorithm does not require the construction of a corresponding Voronoï diagram (usually by employing the Quick Hull algorithm) which typically is a highly computationally demanding task, especially when performed in higher dimensions. The algorithm is suitable for both the clipped and periodic variants of Voronoï tessellation and provides a significant speedup in a convenient range of practical usage. For the purposes of approximation of selected scalar properties of Voronoï cells, only the distances of points in sample are evaluated using an adequately fine underlying orthogonal mesh. The algorithm estimates e.g. the volumes and centroids of Voronoï cells, radii and coordinates of centers of the corresponding Delaunay hyper-circles. The paper also provides quite accurate explicit error estimators for the extracted attributes due to rasterization and thus the user is given a control over the accuracy by selecting an appropriate discretization density. Among numerous fields in which Voronoï diagrams are being utilized, the authors are concerned with optimization of point samples for the Design of Experiments and also with weighting of integration points in Monte Carlo type integration. In these applications, selected scalar topological descriptors of Voronoï diagrams are being repeatedly computed. As the full Voronoï diagram is not of interest but the resulting cell volumes or shape descriptions have to be repeatedly computed, the presented parallel solution seems highly suitable for these applications.

Introduction

The Voronoï diagram has become a classical construct of computational geometry, see [2], [67]. Consider a domain of a general integer dimension, D, containing Ns points (“seeds”, “sites” or “generators”). The task of the Voronoï tessellation is to divide the given space into i=1,,Ns convex polyhedrons of dimension D that encompass regions where all points are closer to ith point than to any other of Ns generating points. A Voronoï tessellation emerges by uniform radial growth from seeds outwards.

Applications of Voronoï diagrams can be found in numerous fields of research and industry (for example biology [42], Thiessen polygons in hydrology [12], [59], chemistry, astrophysics, fluid dynamics [62], computational physics [35], epidemiology [34], medical diagnosis [57], engineering, geometry [70], informatics [46], etc.) The geometric sense of Voronoï diagrams is used for solving tasks like searching for the nearest neighbor, the largest empty circle or estimating the roundness of an object. Technical applications include e.g. construction of meshes in bordered or periodically repeated domains, see [73] and [72].

This paper is motivated by applications of Voronoï tesselation in the field of Design of Experiments. There, the design domain (usually transformed into a convex unit hypercube [0, 1]D) is being filled by a set of Ns distinct points (a “sample” or “design”). The selection of optimal sampling points from a unit hypercube is an old problem that finds many applications in science and industry [5], [11], [18], [19], [21], [24], [25], [29], [36], [40], [43], [50]. A Voronoï diagram can be used for evaluation of uniformity of a point sample and possibly for further optimization of the point layout, see [33], [52], [53]. Additionally, it has been also proposed to use the volume of each cell in the Voronoï diagram as the weight of the corresponding integration point as used in the weighted numerical integration of Monte Carlo type [68], [69]. There is a large potential in these applications to utilize an accelerated extraction of selected scalar descriptors of Voronoï diagrams.

The construction of the exact Voronoï diagram is a computationally demanding task. In the worst cases, the solution complexity approaches O(Ns2), see e.g. the QuickHull algorithm [4]. Moreover, the computational time and the storage demands grow exponentially with the dimension of the design domain, D. On top of these computations, the eventual enumeration of volumes of all D-dimensional cells, see [10], has to be executed.

For many applications, the exact geometrical description of Voronoï cells is not needed; only selected attributes of the regions are of interest. The authors of [48] proposed to solve the largest empty sphere problem (a byproduct of Voronoï tessellation) in the multidimensional space by using a popular stochastic search approach, evolutionary algorithm (EA). They followed the path presented in [39]. In the present paper, a different approach to approximation is selected. The presented contribution proposes a fast parallelized implementation of approximation of volumes of the Voronoï cells without the need of constructing the Voronoï diagram itself. This approach turns out to be useful in cases where the actual shape of the diagram is not of interest but its certain scalar properties have to be repeatedly computed for optimization purposes.

The presented contribution proposes to conduct a rasterization of the hyper-dimensional design space using an underlying mesh of nodes. Because a mesh of equal and independent nodes is highly suitable for a kind of massively parallel execution, the Nvidia CUDA [51] platform is utilized for implementation of the solution.

Section snippets

Discretization approach

The essence of the solution method is based on division of the entire design space into evenly distributed sub-regions. Further used is the division using a regular orthogonal grid of cells. The side of a single cell of such a hyper-dimensional orthogonal grid is considered as 1/ng, where ng is number of grid cells along each axis of the hypercubical domain. The total number of grid cells Ng is then simply Ng=(ng)D.

Cells of a regular orthogonal grid exhibit advantageous properties of equal

Attributes extracted from the rasterized Voronoï diagram

As soon as the index of the closest site is obtained for each of the grid nodes, it is possible to extract certain desired attributes of such a rasterized Voronoï diagram. For purposes of the use case of sample optimization employed by the authors, attributes discussed below are approximated.

Analysis of approximation errors

The proposed approximation of Voronoï tessellation by “rasterization” or “pixelization” of the domain introduces errors in the approximated properties. It is important to know the loss of information and accuracy inherent in the grid, and its relationship to the grid cell size. The grid cells form hypercubes of edge length a=1/ng and their volume is v=aD=1/ngD. The maximum distance within a single grid cell (sometimes referred to as the “pixel”) is its spatial diagonal:ld=aD=Dng.

Implementation of parallel solution

The actual parallel implementation consists of several steps for which the approach to parallel execution differs to some extent. Generally, however, solution of all tasks described in Section 3 does exhibit a rather insignificant arithmetic workload. The dominant scale of the problem is clearly the grid of nodes. The number of nodes in grid rises steeply with (i) the side of grid, ng, that is considered as ng=3Ns, and more so with (ii) the dimension of the domain, D.

With the size of grid being

Performance and speedup

The execution performance of the presented parallel approximation is compared to the standard QuickHull algorithm [4]. The execution times are compared in Fig. 7 left. The solution speedup attained by the parallel solution is shown in Fig. 7 right.

The results are displayed in dependence on the number of sites, Ns, and dimensions, D. In case of the parallel approximation, the side of grid is set as ng=3Ns. Therefore the grid contains Ng=(3Ns)D cells.

In an average case, the complexity of the

Conclusion

The presented paper introduced a parallelized implementation of a numerical approximation of selected scalar properties of a hyper-dimensional Voronoï  diagram. The developed approximation algorithm is motivated by the field of interest of the authors that is optimization of uniformity of point samples for Design of Experiments.

The notion of parallel processing of an underlying raster of nodes crucially relies on the use case of the algorithm. In the case of Latin Hypercube (LH) point samples,

Acknowledgement

This work has been supported by the Czech Science Foundation under project No. GA16-22230S. Additionally the work has also been supported by Specific university research project MSMT No. FAST-J-18-5254 (The Ministry of Education, Youth and Sports of Czech Republic). The authors thank Dr. Václav Sadílek for performing the analysis of the average number of facets of Poisson–Voronoï tessellation.

References (73)

  • P.R. Lloyd

    Quantisation error in area measurement

    Cartogr J

    (1976)
  • R. Miller et al.

    Geometric algorithms for digitized pictures on a mesh-connected computer

    IEEE Trans Pattern Anal Mach Intell

    (1985)
  • L. Pronzato

    Minimax and maximin space-filling designs: some properties and methods for construction

    Journal de la Société Française de Statistique

    (2017)
  • D. Sánchez-Gutiérrez et al.

    Fundamental physical cellular constraints drive self-organization of tissues

    EMBO J

    (2016)
  • M. Vořechovský et al.

    Application of Voronoï weights in Monte Carlo integration with a given sampling plan

  • T. Asano et al.

    Constant-work-space algorithms for geometric problems

    J Comput Geom

    (2011)
  • F. Aurenhammer

    Voronoi diagrams—a survey of a fundamental geometric data structure

    ACM Comput Surv (CSUR)

    (1991)
  • M. Avellaneda et al.

    Weighted monte carlo: a new technique for calibrating asset-pricing models

    Int J Theor Appl Finance

    (2001)
  • C.B. Barber et al.

    The quickhull algorithm for convex hulls

    ACM Transactions on Mathematical Software (TOMS)

    (1996)
  • D.R. Bellhouse

    Area estimation by point-counting techniques

    Biometrics

    (1981)
  • A.K. Bregt et al.

    Determination of rasterizing error: a case study with the soil map of the netherlands

    Geogr Inf Syst

    (1991)
  • W. Brostow et al.

    Voronoi polyhedra and delaunay simplexes in the structural analysis of molecular-dynamics-simulated materials

    Phys Rev B

    (1998)
  • C. Buchta et al.

    Stochastical approximation of convex bodies

    Mathematische Annalen

    (1985)
  • B. Büeler et al.

    Exact volume computation for polytopes: a practical study

    Polytopes–combinatorics and computation

    (2000)
  • P. Burrough

    Principles of geographical information systems for land resources assessment

    Monographs on Soil and Resources Survey (Book 12)

    (1986)
  • C.F. Carver et al.

    Vector to raster conversion error and feature complexity: an empirical study using simulated data

    Geog Inf Syst

    (1994)
  • T.M. Chan

    Semi-online maintenance of geometric optima and measures

    SIAM J Comput

    (2003)
  • Liao Yang et al.

    An equal area conversion model for rasterization of vector polygons

    SCIENCE CHINA Earth Sci

    (2007)
  • L.P. Chew et al.

    Voronoi diagrams based on convex distance functions

    Proceedings of the first annual symposium on Computational geometry - SCG ’85, held in Baltimore, Maryland, USA

    (1985)
  • C. Fabri A. Dehne et al.

    Scalable and architecture independent parallel geometric algorithms with high probability optimal time

    Proceedings of 1994 6th IEEE Symposium on Parallel and Distributed Processing, held in Dallas, TX, USA (26–29 Oct. 1994)

    (1994)
  • Faber. Vance et al.

    Centroidal Voronoi tessellations: Applications and algorithms

    SIAM Rev

    (1999)
  • L. Filipe et al.

    Efficient incremental sensor network deployment algorithm

    In Brazilian Symposium on Computer Networks

    (2004)
  • Y.S. Frolov et al.

    The accuracy of area measurement by point counting techniques

    Cartograph J

    (1969)
  • A.E. Galashev et al.

    Computer simulation of the absorption of CO2 molecules by water cluster: 2. the microstructure

    Colloid J

    (2005)
  • M.F. Goodchild et al.

    Estimation from grid data: the map as a stochastic process

  • M.F. Goodchild

    Fractals and the accuracy of geographical measures

    Math Geosci

    (1980)
  • Cited by (8)

    • An incremental algorithm for simultaneous construction of 2D Voronoi diagram and Delaunay triangulation based on a face-based data structure

      2022, Advances in Engineering Software
      Citation Excerpt :

      For example, Guibas, Knuth, and Shariri [8] Suggested an algorithm with O(logn) overall complexity to find the triangle containing the newly added point by storing the old structures of the DT and connecting the old and new triangles by pointers from the old triangles to the new triangles. With the rapid development of computer hardware, parallel computing based on DT has attracted more attention in recent years [47–50]. But the need for multi-core CPUs is the limitation of the widespread use of the DT or the construction of VD in parallel [40].

    • First-principles study of the diffusion of Li in bcc Fe

      2019, Fusion Engineering and Design
      Citation Excerpt :

      When a Li atom is in the <111>Fe-Li interstitial site, it is under the highest compressive load in <111> direction. To assess the hydrostatic pressure around the Li atom, we calculated the volume of the Voronoi tessellation [35] around the Li atom, where every point inside the cell is closer to the Li atom rather than other Fe atoms. This volume of the Voronoi tessellation is not the exact size of the Li atom but reflects the space around the Li atom and has been widely used in evaluating atomic characteristics, e.g. atomic stress [36].

    View all citing articles on Scopus
    View full text