A Novel Point Inclusion Test for Convex Polygons Based on Voronoi Tessellations
AA Novel Point Inclusion Test for Convex Polygons Based on VoronoiTessellations
Rahman Salim Zengin a,1, ∗ , Volkan Sezer a,2 a Mechatronics Education & Research Center (MEAM), Istanbul Technical University Maslak Campus, Istanbul, Turkey
Abstract
The point inclusion tests for polygons, in other words the point in polygon (PIP) algorithms are fundamentaltools for many scientific fields related to computational geometry and they have been studied for a long time.PIP algorithms get direct or indirect geometric definition of a polygonal entity and validate its containmentof a given point. The PIP algorithms which are working directly on the geometric entities derive linearboundary definitions for the edges of the polygon. Moreover, almost all direct test methods rely on the twopoint form of the line equation to partition the space into half-spaces. Voronoi tessellations use an alternateapproach for half-space partitioning. Instead of line equation, distance comparison between generator pointsis used to accomplish the same task. Voronoi tessellations consist of convex polygons which are definedbetween generator points. Therefore, Voronoi tessellations have become an inspiration for us to develop anew approach of PIP testing specialized for convex polygons. Essential equations to the conversion of aconvex polygon to a voronoi polygon are derived along this paper. As a reference, a very standard convexPIP testing algorithm, the sign of offset , is selected for comparison. For generalization of the comparisons the ray crossing algorithm is used as another reference. All algorithms are implemented as vector andmatrix operations without any branching. This enabled us to benefit from the CPU optimizations of theunderlying linear algebra libraries. All algorithms are tested for three different polygon sizes and varyingpoint batch sizes. Overall, our proposed algorithm has performed better with varying margin between 10%to 23% compared to the reference methods.
Keywords: point inclusion test, point in polygon, convex polygon, voronoi tessellations
1. Introduction
Various point inclusion tests [1] are used in manyapplications[2], including planning for autonomousdriving [3], geographical information systems [4, 5,6], and computer graphics [7, 8]. Any improvementon the efficiency of point inclusion tests will providea direct benefit to the mentioned areas.When autonomous driving related planning ap-plications are considered, planning is mostly donein a 2D space. Collected real time sensor data, es-pecially Lidar based point cloud data is mapped tothe 2D space where the planning operation is per-formed. Collision check is one of the most critical ∗ Corresponding author
Email addresses: [email protected] (RahmanSalim Zengin ), [email protected] (Volkan Sezer ) ORCID: https://orcid.org/0000-0002-3104-4677 ORCID: https://orcid.org/0000-0001-9658-2153 components of motion planning. Several simplifica-tions on collected data and vehicle representation isrequired to make it efficient. Modeling the vehicleas a circle or combination of several circles is oneof the widely used techniques for collision check.Although this simplification works well for most ofthe situations, there is always an accuracy problemdepending on the number of circles that are used tomodel the vehicle [3].In order to make a more accurate collision check,vertical footprint of the car can be modeled as aconvex polygon in a suitable manner by simply find-ing a convex hull bounding the car. In order tomake a real time motion planning, efficient collisioncheck algorithms that are able to test big batches ofpoints against the convex polygon model of the carare needed. Even though there are simple and wellknown algorithms, we propose an alternative algo-
Preprint submitted to Applied Mathematics and Computation August 11, 2020 a r X i v : . [ c s . C G ] A ug ithm based on Voronoi approach to accomplish thesame task.Geographical information systems [4, 5, 6] are an-other field that relies on point inclusion tests. Itis used to process large databases of cartographicdata. Measurements taken in the field must bematched with the prior information related to thearea. For that reason, point inclusion tests are runon big databases using the measurements.Another field in which the point in polygonqueries are being actively used is, computer graph-ics [7, 8]. A scene is composed of many many ob-ject models which are composed of polygons. Forvisualization on the screen, proper rasterization ofthe geometric data is needed. To match the pixelson the screen with the geometric data according tothe defined camera model, polygons are mappedto the screen plane. Then membership of everyscreen pixel is determined via point inclusion test-ing, so that pixels can be painted properly. Geo-metric model of the object which is subject to pointinclusion testing is mostly known priorly and notchanging. Considering this fact, instead of usingpolygonal model of the object directly, a prepro-cessed and simpler to use equivalent definition ofthe object can be utilized.Voronoi tessellations are simple, understandableand well known. They consist of convex polygonsand there is a huge literature related to voronoi tes-sellations. Simplest point inclusion tests are basedon line equations and point to line distance calcula-tions. Conversion of a convex polygon to a voronoipolygon has the advantage of using only point topoint distance calculations. Required equations forthe conversion of a convex polygon to a voronoipolygon are derived step by step throughout thispaper. These derivations are directly used for theimplementation of the proposed algorithm.For completeness, we summarize two simple andwell known point inclusion approaches, and thenused these as a comparison to our proposed method.In order to compare the algorithms in an equal man-ner, all algorithms are implemented, using vectorand matrix operations instead of simple loops, withthe help of the related libraries. In this way, com-putations are handled more efficiently.For the relatively smaller size of point batchessuch as 1024 points, due to the conversion over-head our proposed algorithm gives the worst re-sults. However, as processed batch size increasesour algorithm takes the lead. Processing time ofthe new algorithm becomes 23% better for 5 edges, 13% better for 8 edges and 11% better for 11 edges.The structure of the paper is as follows: In (Sec-tion 2) two reference algorithms are mentioned andthe notation which is used throughout the paper isgiven. In (Section 3), conversion of a convex poly-gon to a voronoi polygon is described and requiredequations are derived. In (Section 4), the pointinclusion testing procedure using the generators isdescribed. In (Section 5), expected performance ofour proposed algorithm is discussed. In (Section6.1), for a certain generated test data our proposedalgorithm is compared against the sign of offset al-gorithm to prove its correctness. In (Section 6.2),experimental setup is described, experimental re-sults are shared and discussed.
2. Background
The use of line equations for linear boundariesare explained in this section via two standard pointinclusion algorithms. method [9, 1, 10] is the goldenstandard of the point inclusion tests for polygons.It can be used for point inclusion testing of simplepolygons. As can be seen on the (Figure 1) a raydirected to the + x direction is used to count cross-ings of the ray and the polygon. If the ray crossesthe polygon edges in odd numbers, it is inside, oth-erwise it is outside.All edges of the polygon are traversed andchecked whether they are on the same y level ofthe point. If applicable, line equation in the twopoint form [11] is used to determine the half-planethe point is present. For a + x going ray it must beon the left half-plane of the line. If so, it is countedas a crossing.The pseudocode of the branchless ray crossingimplementation which is used for experiments isgiven in (Algorithm 1). [5, 1] method is the simplestpoint in polygon algorithm specialized for convexpolygons. A point in a convex polygon given in(Figure 2). For subsequent vertices of the polygon,the offset of the point from the line passing throughthe edge is calculated using the two point form ofline equation. If the offset of the point has thesame sign for all edges of the polygon, the point2 igure 1: Ray crossing method Function
CrossingInclusion
Data: // Q i : Vertices, (2 × n ) Q i ←− (cid:2) q · · · q n (cid:3) // X : Points, (2 × m ) X ←− (cid:2) x · · · x m (cid:3) Result:
IsIn : Boolean, ( m, ) begin // Q j : Rolled vertices, (2 × n ) Q j ←− (cid:2) q n , q · · · q n − (cid:3) // Edges in Y range, ( n × m ) InRange ←− ( Q i > X ) ⊕ ( Q j > X ) // LHS & RHS of the line equation, ( n × m ) LHS ←− ( X − Q j ) ◦ ( Q i − Q j ) RHS ←− ( X − Q j ) ◦ ( Q i − Q j ) // Is edge going up, ( n, ) U p ←− Q i > Q j // Is point on the left, ( n × m ) OnLef t ←− U p ? (
LHS > RHS ) :(
RHS > LHS ) Crossing ←− InRange ∧ OnLef tIsIn ←− M od ( (cid:80) i Crossing ij ) = 0 endAlgorithm 1: Ray crossing point inclusion test
Figure 2: Sign of offset method is inside. Otherwise, if the two subsequent edgesgive different signs for the offset of the point, it isoutside.The pseudocode of the implemented branchlessvariant of the algorithm is given in (Algorithm 2).
Function
SignOfOffsetInclusion
Data: // Q i : Vertices, (2 × n ) Q i ←− (cid:2) q · · · q n (cid:3) // X : Points, (2 × m ) X ←− (cid:2) x · · · x m (cid:3) Result:
IsIn : Boolean, ( m, ) begin // Q j : Rolled vertices, (2 × n ) Q j ←− (cid:2) q n , q · · · q n − (cid:3) // LHS & RHS of the line equation, ( n × m ) LHS ←− ( X − Q j ) ◦ ( Q i − Q j ) RHS ←− ( X − Q j ) ◦ ( Q i − Q j ) // Sign test, ( n × m ) D ←− LHS < RHS // Are all same sign
IsIn ←− M od n ( (cid:80) i D ij ) = 0 endAlgorithm 2: Sign of offset point inclusion test
For simplicity and clearance, definitions relatedto voronoi tessellations [12] are slightly modifiedand adapted.Throughout this paper, only 2-dimensional eu-clidean space, R is considered. Boldface denotesa vector, such as x = ( x , x ) T . Superscript T n vertices, vertices of the polygon are de-noted with additional indexes, such as q i , q j , where i, j = { , . . . , n } and i (cid:54) = j where edges consid-ered. The set of vertices of the voronoi polygon is Q = { q , . . . , q n } .A voronoi polygon is a convex region defined byan inner generator point and some outer generatorpoints such as, V ( p ) = { x | (cid:107) x − p (cid:107) ≤ (cid:107) x − p k (cid:107)∀ k ∈ { , . . . , n }} (1)where V ( p ) denotes voronoi polygon related to thegenerator point p .A generator point p k belongs to the set of gen-erator points P of the voronoi polygon. Inner gen-erator point is always indexed as p independentof the polygonal edge count n . For every edgeof the voronoi polygon there is an outer gener-ator point so that the set of generator points is P = { p , p , . . . , p n ).Edges are equidistant set of points between innergenerator and outer generators. Precisely, e k = { x | (cid:107) x − p (cid:107) = (cid:107) x − p k (cid:107)} (2)where k ∈ { , . . . , n } . The set of edges of the V ( p )can be denoted as E = { e , . . . , e n } .The whole set of edges of the voronoi polygon iscalled boundary and it is denoted related to theinner generator point of the region encircled as ∂V ( p ). Although a voronoi graph has multiplepolygonal regions, throughout this study, we areonly interested in defining a single voronoi polygon.For two subsequent vertices q i , q j of a polygon,two point form of the line equation [11] can be writ-ten as( x − q i )( q j − q i ) = ( x − q i )( q j − q i ) (3)where the point, x = ( x , x ) T .
3. Conversion of convex polygons to voronoipolygons
Vertices of a convex polygon ( q i on Figure 3) canbe taken as the vertices of a voronoi polygon, andedges of a convex polygon ( e k on Figure 3) can betaken as the boundary of a voronoi polygon.Because determination of generator points ( p k on Figure 3) is only constrained by ∂V ( p ), it isfree to choose any internal point as p . But to Figure 3: Generators, vertices and edges of a voronoi polygon distribute outer generators homogeneously and topreserve symmetry, when the polygon is symmet-ric, and to have a guaranteed point inside, the cen-troid of the polygon is used as the inner point. Thecentroid of a convex object always lies inside, be-cause any line passing through the centroid can onlycross the edges of the convex polygon at two points.Then, outer generator points can be found accord-ingly.As can be seen on (Figure 3) placement of genera-tor points determines not only ∂V ( p ) but also theedges going to the infinity between outer generatorpoints. However, our problem is only constrainedon ∂V ( p ).Having n edges or lines passing through n cou-ples of vertices constrained on ( n + 1) generatorpoints, gives us freedom of choosing one of the gen-erator points. Although setting any of the genera-tor points sets all the others, setting the inner gen-erator is more reasonable, because all the edges aredefined depending upon it.The centroid of a polygon [13] can be calculatedas follows:Let Q be a cyclically ordered set of polygon ver-tices and q i , q j are subsequent vertices accordingly.Summation over Q , A = 12 (cid:88) Q det[ q i q j ] (4) p = 16 A (cid:88) Q ( q i + q j ) det[ q i q j ] (5)4 p k q j q � e k ax + bx + c = 0bx - ax + d = 0 x Figure 4: Finding outer generator of an edge gives the area (4) and centroid (5) of the polygon.The pseudocode of the centroid calculation is givenin (Algorithm 3).
Function
CalculateCentroid
Data: // Q i : Vertices, (2 × n ) Q i ←− (cid:2) q · · · q n (cid:3) Result: µ : Centroid, (2 , ) begin // Q j : Rolled vertices, (2 × n ) Q j ←− (cid:2) q n , q · · · q n − (cid:3) // A : Partial areas, ( n, ) A ←− Q j ◦ Q i − Q i ◦ Q j a ←− (cid:80) A // a : Area, (1 , ) µ ←− (cid:0) ( Q i + Q j ) A (cid:1) / (6 a ) endAlgorithm 3: Calculation of the centroidStandard form equation of the line passingthrough an edge can be derived from two point formequation. q i and q j are two vertices of the edge e k . x is a point of the edge. As defined in (2) p and p k are two points, equidistant to the e k . The linepassing through p and p k is perpendicular to (3).Solving x for two equations gives x = (cid:18)(cid:20) b k − a k b k − a k b k a k (cid:21) p − c k (cid:20) a k b k (cid:21)(cid:19) a k + b k (6) where a k = q i − q j b k = q j − q i c k = − ( a k q i + b k q i ) p and p k are equidistant to x . Writing thisequation and leaving p k alone on the left hand sidegives p k as p k − x = x − p ⇒ p k = 2 x − p (7)By substituting x into (7), outer generator pointscan be found as p k = (cid:18)(cid:20) b k − a k − a k b k − a k b k a k − b k (cid:21) p − c k (cid:20) a k b k (cid:21)(cid:19) a k + b k (8)The generator calculation procedure is given in(Algorithm 4). Function
CalculateGenerators
Data: // Q i : Vertices, (2 × n ) Q i ←− (cid:2) q · · · q n (cid:3) Result: P : Generators, (2 × ( n + 1)) begin P , ←− CalculateCentroid ( Q i ) // Q j : Rolled vertices, (2 × n ) Q j ←− (cid:2) q n , q · · · q n − (cid:3) a ←− Q i − Q j // ( n, ) b ←− Q j − Q i // ( n, ) c = − ( a ◦ Q i + b ◦ Q i ) // ( n, ) // W : (2 × × n ) W ←− (cid:20) b − a − ab − ab a − b (cid:21) d ←− (cid:80) j W ijk P j // (2 × n ) e ←− − c ◦ (cid:20) ab (cid:21) // (2 × n ) P , n +1) ←− ( d + e ) (cid:11) ( a + b ) endAlgorithm 4: Calculation of generators
4. Point inclusion test via generator points
After the set of generators P has been found, thepoint inclusion test is simply testing the conditionprovided in (1).5 unction VoronoiInclusion
Data: // Q : Vertices, (2 × n ) Q ←− (cid:2) q · · · q n (cid:3) // X : Points, (2 × m ) X ←− (cid:2) x · · · x m (cid:3) Result:
IsIn : Boolean, ( m, ) begin // P : Generators, (2 × ( n + 1)) P ←− CalculateGenerators ( Q ) // ∆ : Differences, (2 × ( n + 1) × m )∆ ←− X − P // M : Metrics, (( n + 1) × m ) M jk ←− (cid:80) i ∆ ijk ∆ ijk IsIn ←− M ≤ M k , ∀ k ∈ { , . . . , ( n + 1) } endAlgorithm 5: Voronoi point inclusion testThe ordinary distance metric for the definition ofthe voronoi polygon is euclidean distance or equiv-alently
L2 norm . To test the inclusion of a randompoint, its distances to all generators are calculated.If it is closest to the generator p it is inside of thepolygon. Otherwise it is outside of the polygon.Ordinarily, calculating the L2 norm of a vector(9) takes squaring, summing and then square root-ing of the vector components. (cid:107) x (cid:107) = (cid:113) x + x (9)However squaring of both sides of (1) does notchange the order of distances because squaring is amonotonic operation. V ( p ) = { x | (cid:107) x − p (cid:107) ≤ (cid:107) x − p k (cid:107) } (10)The square root and the square vanish, whenthese are applied together. Then equation (10) be-comes V ( p ) = { x | ( x − p ) T ( x − p ) ≤ ( x − p k ) T ( x − p k ) } (11)The derived simplification (11) is an alternateway of distance comparison. It improves the per-formance of comparisons and preserves the order ofdistances.The pseudocode of the proposed point inclusiontesting algorithm is given in (Algorithm 5)
5. Algorithm analysis
The calculation of polygon centroid takes O ( n )time when it is done sequentially. Similarly, outergenerator point calculations have time complexityof O ( n ). But considering Single Instruction Multi-ple Data (SIMD) capabilities of modern CPUs, forsmall sizes of n computations will be optimized tobe done with time complexity of O (1).For n vertices and m points; ( n + 1) m distancecalculations are done. Then using distances to theinner centroid as a reference, the number of dis-tance comparisons to be made is nm . Conversionrelated computations are done initially and are in-dependent of the batch size of processed points. Asthe batch size m of processed points increases, theconversion cost becomes less effective on the overallcomputational cost.In practice, for determination of the status of apoint, doing all computations and comparisons isnot always needed. If the point under test is foundto be closer to an outer generator, this breaks the ∀ condition of (1). An early break opportunity ariseshere for a sequential implementation of the algo-rithm.
6. Experimental results and discussion
To test correctness of the proposed point inclu-sion algorithm, random test points are sampled(Figure 5) around the polygon. The set of gen-erators for the tested polygon are also plotted.Inclusion test results of the sign of offset algo-rithm are used as the known ground truth refer-ence. For the same test set, both algorithms gavethe same results. The correctness of the proposedalgorithm is proved via this testing procedure. Thecorrectness of the algorithm can be seen in (Figure5) by looking at different coloring of the dots insideand outside.
In order to make a fair comparison, calculationsare performed for all vertices, edges or generatorsetc. every time for each relevant method. Thus,experimental results reflect theoretical complexityanalyses well.The CPU used for the experimentation is In-tel(R) Core(TM) i7-7700 running at 3.60GHz fre-quency. The system has 32G of RAM.6
Figure 5: Correctness test of the proposed algorithm
For ease of reproducibility, all implementationsare done using Python [14] and related libraries [15,16]. The source code [17] to reproduce the resultsis shared.(Table 1) provides information about the process-ing time of each method for various types of con-vex polygons, using different number of test points.Per point processing times are calculated via divid-ing batch processing time into number of points inthat batch. Best methods for each test scenario areshown bold. The difference column in the table iscalculated between the best two results in a rowaccording to the formula (12) below:
Dif f erence % = 2 ndBest − stBest ndBest ×
100 (12)Ray crossing is a more general algorithm whichcan test the inclusion of points for simple polygons.It is provided as one of the references here since itis the golden standard of point inclusion tests forbenchmarking the point in polygon algorithms. Asit is a more general and complicated approach, itis slower than the algorithms specialized for convexpolygons.On the other hand the sign of offset algorithm,which is specialized for convex polygons, gives bestresults for small problem sizes. Conversely, the pro-posed Voronoi based method is the slowest one forrelatively smaller problem sizes since it has a con-stant conversion cost. On the other hand, with in- P e r p o i n t p r o c e ss i n g t i m e ( n s ) Tests for 5 edges ray crossingsign of offsetvoronoi10 P e r p o i n t p r o c e ss i n g t i m e ( n s ) Tests for 8 edges ray crossingsign of offsetvoronoi0 200000 400000 600000 800000 100000010 P e r p o i n t p r o c e ss i n g t i m e ( n s ) Tests for 11 edges ray crossingsign of offsetvoronoi
Figure 6: Tests for 5, 8, and 11 edges able 1: Per point processing times (ns) AlgorithmsRay Sign of DiffEdges Points Crossing Offset Voronoi %1024 161.49 ray crossing following our algorithmand lastly the sign of offset algorithm because ithas least temporary storage requirements.
7. Conclusion
A systematic approach to convert a convex poly-gon to a voronoi polygon is developed throughoutthis work. As a meaningful internal generator pointselection scheme, centroid calculation of a polygonis chosen. The equations related to the centroidcalculation are given consecutively. After that, theequations required to calculate outer generators inrelation to the inner generator and the vertices ofthe convex polygon, are derived.In order to demonstrate the advantages of ourproposed algorithm, it is implemented as only vec-tor and matrix operations without including anyprogram branching. Reference algorithms are alsoimplemented in a similar way. Certain tests are car-ried out to show that our proposed algorithm notonly works properly but also it has a clear perfor-mance advantage over the reference algorithms forlarge size of point batches.Conversion of a convex polygon to a voronoi poly-gon takes constant time according to the pointbatch size. It depends only on the number of edgesof the polygon. If the geometry is known to beconstant prior to the use, voronoi equivalent of theconvex polygon can be calculated in advance. Bothpolygon vertices and generator points can be storedtogether in a database with only about 2 × increaseof the storage capacity needed.Testing vertices of a convex polygon against gen-erators of another one, collision detection betweenconvex polygons can be accomplished. Performingthis testing in a mutually inclusive way between twoconvex polygons reduces errors.The purpose of this paper is showing applicabil-ity of generator points based boundary definitionand determination to the point inclusion testing.Besides, a novel algorithm is matched with priorknowledge in the field. Thus, the methods devel-oped here can be extended to nonconvex polygons.The generator based boundary checks can be ap-plied to prior point in polygon algorithms. Acknowledgements
This work was supported by the TurkishScientific and Technological Research Council(TUBITAK) under project no. 118E809.8 eferences [1] E. Haines, Point in Polygon Strategies, in:Graphics Gems, Elsevier, 1994, pp. 24–46. doi:10.1016/B978-0-12-336156-1.50013-6 .URL https://linkinghub.elsevier.com/retrieve/pii/B9780123361561500136 [2] D. Alciatore, R. Miranda, A winding number and point-in-polygon algorithm, Glaxo Virtual Anatomy ProjectResearch Report, Department of Mechanical Engineer-ing, Colorado State University (1995).[3] J. Ziegler, C. Stiller, Fast collision checking for intelli-gent vehicle motion planning, in: 2010 IEEE IntelligentVehicles Symposium, 2010, pp. 518–522, iSSN: 1931-0587. doi:10.1109/IVS.2010.5547976 .[4] T. K. Peucker, N. Chrisman, Cartographic DataStructures, The American Cartographer 2 (1) (1975)55–69. doi:10.1559/152304075784447289 .URL [5] S. Nordbeck, B. Rystedt, Computer cartography point-in-polygon programs, BIT Numerical Mathematics7 (1) (1967) 39–64. doi:10.1007/BF01934125 .URL https://doi.org/10.1007/BF01934125 [6] P. A. Longley, M. F. Goodchild, D. J. Maguire, D. W.Rhind (Eds.), Geographical Information Systems: Prin-ciples, Techniques, Management and Applications, 2ndEdition, Wiley, Hoboken, N.J, 2005.[7] E. Haines, FAST RAY CONVEX POLYHE-DRON INTERSECTION, Graphics Gems II(1991) 247–250Publisher: Morgan Kaufmann. doi:10.1016/B978-0-08-050754-5.50053-0 .URL [8] A. S. Glassner, An Introduction to Ray Tracing, Else-vier, 1989, google-Books-ID: t3XNCgAAQBAJ.[9] M. Shimrat, Algorithm 112: Position of point relativeto polygon, Communications of the ACM 5 (8) (1962)434. doi:10.1145/368637.368653 .URL http://portal.acm.org/citation.cfm?doid=368637.368653 [10] W. R. Franklin, PNPOLY - Point Inclusion in PolygonTest - WR Franklin (WRF) (2020).URL https://wrf.ecse.rpi.edu//Research/Short_Notes/pnpoly.html [11] E. W. Weisstein, ”Two-Point Form.” From MathWorld– A Wolfram Web Resource, library Catalog: math-world.wolfram.com Publisher: Wolfram Research, Inc.URL https://mathworld.wolfram.com/Two-PointForm.html [12] A. Okabe, Spatial tessellations: concepts and applica-tions of Voronoi diagrams, 2nd Edition, Wiley seriesin probability and statistics, Wiley, Chichester ; NewYork, 2000.[13] G. Bashein, P. R. Detmer, Centroid of a Poly-gon, in: Graphics Gems, Elsevier, 1994, pp. 3–6. doi:10.1016/B978-0-12-336156-1.50010-0 .URL https://linkinghub.elsevier.com/retrieve/pii/B9780123361561500100 [14] G. van Rossum, Python tutorial, technical report CS-R9526, centrum voor wiskunde en informatica (CWI),amsterdam.” (1995).[15] S. van der Walt, S. C. Colbert, G. Varoquaux, TheNumPy Array: A Structure for Efficient NumericalComputation, Computing in Science Engineering 13 (2) (2011) 22–30, conference Name: Computing in ScienceEngineering. doi:10.1109/MCSE.2011.37 .[16] J. D. Hunter, Matplotlib: A 2D Graphics Environment,Computing in Science Engineering 9 (3) (2007) 90–95,conference Name: Computing in Science Engineering. doi:10.1109/MCSE.2007.55 .[17] R. S. Zengin, V. Sezer, volimpro/voronoi pip: REL:v1.0 (Jul. 2020). doi:10.5281/zenodo.3966996 .URL https://doi.org/10.5281/zenodo.3966996https://doi.org/10.5281/zenodo.3966996