 # Overlaid oriented Voronoi diagrams and the 1-Steiner tree problem

Overlaid oriented Voronoi diagrams (OOVDs) are known to provide useful data for the construction of optimal Euclidean 1-Steiner trees. The theoretical time complexity of construction methods exploiting the OOVD is O(n^2), but a computational study has never been performed, and robust constructions for OOVDs have not previously been implemented. In this paper, we outline a numerically stable implementation for constructing OOVDs using tools from the Computational Geometry Algorithms Library (CGAL), and test its performance on random point sets. We then study the effect that the OOVD data has in reducing the complexity of 1-Steiner tree construction when compared to a naive approach. The number of iterations of the main loop of the 1-Steiner algorithm is directly determined by the number of faces in the OOVD, and this appears to be linear for the random inputs we tested. We also discuss methods for processing the OOVD data that lead to a reduction in construction time by roughly a factor of 12.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1. Introduction

The Euclidean Steiner tree problem is to construct the network with shortest total edge length that connects a set of input terminals in the plane. The leading algorithm for this problem, GeoSteiner, has at its core a pruning technique that uses the fact that any Euclidean minimum spanning tree (MST) on a set of points is a subgraph of the relative neighbourhood graph of  [13, 9, 14]. The relative neighbourhood graph is a subgraph of the Delaunay triangulation of , and is therefore related to its dual, the Voronoi diagram of . Optimal Euclidean Steiner trees contain new junctions, or Steiner points, which all have degree 3.

The Euclidean -Steiner problem on is to construct the network with shortest total edge length that connects the points of and uses at most Steiner points. This problem requires a different approach to the original Steiner tree problem for the simple reason that Steiner points of degree can occur in an optimal solution. The success of the GeoSteiner algorithm in efficiently solving large Steiner tree problems is crucially tied (via the Melzak algorithm) to the assumption that all Steiner points are of degree .

The Steiner tree problem is NP-hard, while it is fairly easy to show that the -Steiner tree problem is solvable in polynomial time for constant . That said, it takes significant ingenuity to keep the degree of the polynomial reasonably low. The first result in this direction was a quadratic time algorithm for the case described by Georgakopoulos and Papadimitriou in 1987 . In that paper they demonstrated the utility of another kind of Voronoi diagram, the overlaid oriented Voronoi diagram.

Recall that the classical Voronoi diagram of a set essentially records, for each point in the plane, the closest point in . An oriented Voronoi diagram (OVD) for and a given cone records, for each point in the plane, the closest point of within the translate of based at . We think of the cone informally as a ‘range of directios’, e.g. . The overlaid oriented Voronoi diagram (OOVD) is the overlay of six OVDs corresponding to six equiangular cones that partition the plane. Thus it records, for each point , the closest point of in each of six different ranges of directions covering the whole plane. Examples of OVDs and an OOVD can be seen in Figure 1.

Monma and Suri  independently discovered the power of the OVD when studying so-called ‘transitions’ in minimum spanning trees. The question they asked is: how is a minimum spanning tree topology affected by perturbations of the the input terminals? The OVD structure allowed them to bound the number of minimum spanning tree topologies that can result when one of the terminals is perturbed. The problem of efficiently predicting topological changes of minimum spanning trees is important for sensitivity analysis, for instance in cases where the input terminal locations are not known with certainty. Models for transitions in MSTs also capture the topological changes in dynamic networks where, for example, input terminals are mobile.

A number of papers have since provided further theoretical results about the algorithmic construction of OVDs and their application to network design. In  an sweep algorithm is presented for constructing the OVD. Note however, this does not provide an overall reduction in complexity of OOVD construction, since the overlay process still requires time. In  OOVDs were applied to the more general problem of optimal -Steiner tree construction in various planar norms and for various cost functions on the edges, yielding an algorithm of complexity , for some function . The OOVD structure has also been applied to the design of survivable networks .

Despite the importance of OOVDs, we are not aware of any practical implementations or experimental studies of OOVD construction. In Section 2 we describe the mathematical construction of the OVD that we use, discuss the CGAL code that implements it and creates the OOVD, and report on the performance of the code on random input sets. We find that the code is robust and appears to run in much better than the worst case quadratic time complexity on these inputs. The complexity of the OOVD itself seems essentially linear.

The time algorithm for the 1-Steiner problem presented in  is theoretical and hasn’t been refined or implemented practically. In Section 3 we discuss the 1-Steiner problem and the way the OOVD data helps to solve it. We use new pruning methods to reduce the OOVD data that result in roughly a factor of 12 improvement in running time over a naive implementation, and also make some statistical observations about the optimal 1-Steiner trees found during our experiment. The efficiency of the algorithm is closely tied to the complexity of the OOVD, and thus also appears to be much better than the worst case analysis suggests.

## 2. The overlaid oriented Voronoi diagram

Given a set of input points, the oriented Voronoi diagram is a partition of the plane into cells that records, for each point in the plane, the closest input point in a fixed range of directions. To define the OVD formally, let be a fixed cone (based at the origin) and consider a set of input points. Then the cell of the OVD corresponding to the point is

 {x∈R2:x−p∈C & (∀q∈P:x−q∈C⇒|x−p|⩽|x−q|)}

Alternatively, we can define the cell in terms of translated cones based at , and halfplanes of points closer to than . Then the cell corresponding to is

 Cp∖⋃q∈P∖{p}Cq∩Hqp.

Note that unlike Voronoi diagrams, the cells of an OVD are not always convex and do not cover the entire plane.

In the case of the -Steiner problem, we use six OVDs corresponding to six cones of angle that partition the plane. The overlaid oriented Voronoi diagram is the common refinement of these six OVDs. Therefore each cell of the OOVD is the intersection of a cell from each of the six OVDs. Thus it records the closest input point in each of the six ranges of directions determined by the cones.

### 2.1. Constructing the oriented Voronoi diagrams

We use a very general method for calculating Voronoi type diagrams that uses the upper envelope of a set of surfaces in 3 space . Given a set of surfaces in , considered as graphs of functions of and , the upper envelope is a partition of the -plane according to which surface is highest above each point . In other words, it is a planar map with each region corresponding to a subset of the domain for which a particular function has the highest value.

It is well known that the classical Voronoi diagram can be constructed as the upper envelope of a set of planes tangent to the unit paraboloid, [6, 5]. Consider a point . The tangent plane to at has equation . Observe that at another point , the vertical distance between and is

 f(q)−hp(q)=q2x+q2y−2pxqx−2pyqy+p2x+q2x=(qx−px)2+(qy−py)2.

That is, it is the square of the distance from to . Therefore, given a set of points , if one considers the set of tangent planes , then for any point the plane that has the highest value at corresponds to the closest point in to . This means the classical Voronoi diagram of is precisely the upper envelope diagram for the planes .

The OVD can be constructed in the same way by using a cone in each of the tangent planes with apex at . In practice we use a triangle which is this cone truncated at one of the lines bounding the problem domain. After constructing these triangles we calculate the upper envelope to create the OVD. Once the OVDs are constructed, overlaying them is mathematically straightforward.

For reasons that we discuss in Section 3, we also construct the classical Voronoi diagram. We use the same triangle envelope algorithm, but instead of using tangent planes, we create large tangent triangles that cover the entire domain. We will overlay the OOVD and the classical Voronoi diagram to create what we call the refined OOVD.

Figure 1 shows an example of the six OVDs and the OOVD for a set of 10 input points. It also includes the classical Voronoi diagram, the refined OOVD, and the standard OOVD.

### 2.2. Overview of the CGAL code

We decided to use the Computational Geometry Algorithms Library (CGAL) for OOVD construction because it has robust and fast implementations of the main geometric procedures we require, namely an upper envelope function for sets of triangles in and an overlay function that returns the common refinement of two planar maps. These functions are discussed in great detail in .

The following is a brief overview of our CGAL code for OOVD creation. We defer detailed discussion of the code until we are able to release it publicly. From the input points, our CGAL program first creates the triangles required for the upper envelopes. To these triangles we attach data that records the associated input point. Next, the seven envelope diagrams are created. The data from the triangles is passed to the regions of the diagrams, associating each region with the relevant input point. After trimming the diagrams down to the bounding box of the problem domain, we begin overlaying them. The overlay function overlays two diagrams at a time, so we need 6 overlay operations to combine them all. During each overlay operation, the data from the faces of the input diagrams is combined to create the data for each face of the output diagram.

The data is only valid on full dimensional faces, not edges or vertices. The data from the faces of the OOVD, rather than the diagram itself, is what is used in the -Steiner algorithm.

### 2.3. Numerical accuracy and coordinate number type

The validity of the data produced by our program depends on the faces of the OOVD being calculated precisely. Even with rational input points, the OVDs contain a lot of edges with irrational slopes, due to the angles of the cones. Representing the diagrams with rational coordinates (including floating point) resulted in rounding errors that caused extra faces to appear – essentially long thin faces that should have been edges. Since the number of faces directly affects the complexity of the -Steiner algorithm, this was a significant issue.

One solution to the problem is to use coordinates in the quadratic extension field . Geometrically, since computers inevitably use bounded precision to represent reals, one can think of this as using coordinates in a triangular rather than square grid. It can be shown that, for rational input points, all intersections of lines in the OVDs and their overlays can be represented exactly in these coordinates. Fortunately CGAL comes equipped with the ability to use coordinates in a quadratic extension field. Our experiments confirmed that this approach eliminated the errors observed using rational coordinates.

### 2.4. Theoretical complexity of computing the OOVD

Constructing the OOVD is part of the preprocessing, and is independent of for the -Steiner problem. The -Steiner problem can be solved in , while the CGAL algorithms that we use have theoretical running times slightly more than . However this turns out to be a very minor issue, as we will see. For the theoretical bounds on the complexity of the algorithms we use are much lower than the complexity of the -Steiner algorithm.

The geometric complexity (number of faces, edges and vertices) of an OVD is  , however the complexity of constructing the OVD is , using a sweepline type algorithm . In general, the upper envelope based method we use for computing the OVD has time complexity , where is the inverse Ackermann function . This bound is the same as the worst case geometric complexity of the upper envelope of a set of triangles, so there can’t be a faster construction for the upper envelope of triangles in general. The algorithm used in CGAL may have even worse complexity of because it accepts collections of more general surfaces111Actually semi-algebraic surfaces of constant descriptive complexity. . However, as the complexity of an OVD is much less than this bound, we expect that the true running time of the upper envelope algorithm for the special case of OVD construction is significantly better than . Of course it can’t be better than the lower bound in  which is due to a reduction from the sorting problem.

Turning to the task of overlaying the six OVDs, the geometric complexity of the OOVD is in the worst case, and Georgakopoulos and Papadimitriou describe a straightforward algorithm to construct it in time  . They use the observation that we can reduce the question to the problem of overlaying triangulations, since the OVDs can be triangulated with triangles. Even without this simplification, it is known that the complexity of the overlay of a constant number of upper envelope diagrams of surfaces in is  . Thus the CGAL overlay algorithm may again have slightly worse than quadratic complexity in theory.

### 2.5. Experimental setup and results

In practice the CGAL based implementation of OOVD construction worked very effectively.

The OOVD code was tested on uniformly randomly generated point sets in a 10,000 by 10,000 integer grid. The main experiment involved computing OOVDs for sets of size 10, 20, 50, 100, 200 and 500, with 50 instances of each size. In each case we computed the six OVDs, the classical Voronoi diagram, the OOVD, the refined OOVD, and the output was just the data from the faces of the refined OOVD. Figure 2 shows the user time for all these instances and a line connecting the averages. All OOVD computations were performed on a system with an Intel Core i7-8650U CPU at 1.9 GHz (4.2GHz turbo). The system had 16GB of RAM, although very little was used.

The shape of the curve in Figure 2 could suggest that on this kind of example the algorithm runs faster than the worst case time discussed in Section 2.4. However there is also a wide variation in time taken, making it hard to draw strong conclusions. For example, the 500 point examples took between 19 and 25 seconds to compute. Figure 2. User time by problem size (averages of 50 seeds).

Overlaying the classical Voronoi diagram with the OOVD obviously increases the number of faces, and thus the amount of data to be used later. To illustrate the effect of this, we compared the average number of faces in the OOVD with and without the Voronoi overlay222Note however that number of faces in the OOVD alone was only calculated on a smaller set of examples.. The resulting plot in Figure 3 shows a strikingly linear relationship between the number of input terminals and the number of regions in the OOVD. For example the number of cells in the refined OOVD is a little under 50 times the number of terminals. This is perhaps not surprising given the points were selected uniformly at random within the square domain. One could expect that the diagram looks locally similar everywhere regardless of the number of points (except near the domain boundary), so that the typical number of cells near a terminal is fairly constant. Such local behaviour would lead to the linear relationship that was observed.

It is feasible to calculate much larger instances on the same hardware. For example, one run on an input set of 10,000 points was completed in just over 17 minutes. This produced a refined OOVD with 483,000 regions, supporting the apparently linear relation in Figure 3. The computation used almost 3GB of RAM at one point.

In summary, our results thus far suggest that, for the kind of input we used, the number of regions in the OOVD is linear and the time taken for construction is not much worse than linear. Clearly further experiments could be run on larger examples to explore the performance of the code further. Other options for further investigation include using different random distributions of input terminals, using structured input sets (to further test robustness), and finer time analysis to determine which part of the process takes more time, the OVD construction or the overlaying steps.

## 3. 1-Steiner trees as an application of OOVDs

In this section we describe an algorithm, based on the one proposed in , which employs the refined OOVD output in order to construct optimal Euclidean -Steiner trees. We then analyse the improvement in efficiency that this provides over a more naive approach. Figure 4 shows two examples of optimal -Steiner trees. Figure 4. Minimum 1-Steiner trees for (top) 10 terminals, and (bottom) 500 terminals.

Constructing optimal -Steiner trees is a polynomially solvable problem. This follows from two observations. Firstly, an optimal -Steiner tree is an MST on its nodes. Therefore, if the optimal location of the Steiner point is known then an optimal -Steiner tree can be constructed using Prim’s algorithm. Secondly, if the neighbours of a Steiner point are known then the optimal location of the Steiner point, which coincides with the Fermat point of these neighbours, can be found in constant time. This is because a Steiner point never has degree more than in the Euclidean plane , and the Fermat point of or of points can be exactly constructed in constant time .

From the above observations, a naive algorithm for finding an optimal -Steiner tree simply iterates through all sets of three or four terminals in and finds the Fermat point of each set. An MST is then calculated on each set and the shortest resulting tree is selected as the optimal solution. This naive approach has an asymptotic complexity of , assuming that a Euclidean MST can be constructed in time .

Employing the OOVD, the worst case complexity for -Steiner tree construction can be reduced to as follows: first the OOVD is created for the given set of input terminals and an MST on is constructed. Next preprocessing is performed that stores the longest edge on every path in . This data will be used to update , in constant time, to include a Steiner point. For each region of the OOVD and its associated terminals, one selects a subset of at most terminals as neighbours for the Steiner point and then calculates the optimal location of the Steiner point with respect to these neighbours. The MST is then updated in order to contain the Steiner point and its neighbours. Below we describe some of these steps in more detail and also present various refinements and practical improvements to this framework.

### 3.1. Calculating the Fermat point

Suppose that is an an optimal -Steiner tree on with Steiner point , where has neighbour-set for some index set . It is known that is of degree at most four , so we may assume that . Since minimises the sum of distances to its neighbours, is the Fermat point of . Although there is no general algebraic solution for constructing Fermat points, there are simple constructions for three or four neighbours. If then forms a triangle with all internal angles at most and where every pair of edges incident to meet at exactly . When the elements of are in convex position and lies at the intersection of the two diagonals of  .

### 3.2. Processing the OOVD output

Each face of the refined OOVD stores the closest terminal in each of the six cones plus the overall closest terminal. This data is stored as a set

of vectors that have 7 entries which are the indices of the relevant terminals in

(or zero if there is no terminal in that cone), with the last entry for the closest terminal. For a vector we denote the associated face of the refined OOVD by . Thus if a point is a Steiner point in an optimal 1-Steiner tree, the indices of the neighbours of all occur in .

The 1-Steiner algorithm involves iterating through every potential set of neighbours of the Steiner point that is to be added to . We call such a potential neighbour set a bucket. Naively, each subset of or of the first six entries of each needs to be considered as a bucket, yielding up to buckets per face. In this subsection we describe a method of processing the OOVD data to produce a significantly smaller set of buckets .

Firstly, the local geometry of Steiner points allows us to ignore certain combinations of elements in each which cannot correspond to the neighbours of a Steiner point. Take any and suppose is an optimal Steiner point lying in . If is of degree , since each cone has an angle of and the edges incident to a degree- Steiner point meet at , there are only two feasible subsets of size , namely and . If is of degree then it lies at the intersection of the two diagonals of the quadrilateral formed by its neighbours. Therefore the neighbours of are two pairs of terminals that occur within opposite cones and so the feasible subsets of size in are , and . Each face therefore yields at most buckets in .

Moreover, since the Steiner point in an optimal 1-Steiner tree must be adjacent to its closest neighbour in , any bucket derived from can be eliminated if it does not contain . This eliminates one of the 3-entry buckets, and one of the 4-entry buckets, leaving at most buckets per face.

Finally, note that it is possible for the same bucket to arise from multiple faces of the refined OOVD. The final stage of processing consists of removing such repetitions.

### 3.3. Preprocessing for the MST update step

One of the steps in the algorithm for constructing an optimal -Steiner tree consists of updating a MST to include a new point . This step can be done in constant time if suitable preprocessing has been performed. This preprocessing step calculates a minimum spanning tree on , and then for each pair of terminals stores in a table the longest edge on the path between these nodes in the MST. The preprocessing step can be perfomed in a time of [ref GenkStei].

### 3.4. The 1-Steiner algorithm

We now present pseudo-code for the -Steiner algorithm.

The complexity of the -Steiner tree algorithm is . This follows from the fact that Line 1 can be performed in time, since contains buckets. Then, Line 2 requires time and Line requires time. The For loop in Line 4 runs at most times, and in each execution of the loop the Fermat point is constructed and the MST updated in constant time.

### 3.5. Experimental results for 1-Steiner trees

In this section we provide experimental results on the effect that the OOVD face data preprocessing has on the number of buckets that must be considered during -Steiner tree construction. Recall that construction time is directly proportional to the number of buckets considered since updating an MST and constructing a Fermat point are done in constant time. We also provide some statistical observations on properties of the optimal -Steiner trees that arose. Figure 5. Number of buckets by problem size (averages of 50 seeds).

Figure 5 shows the number of buckets generated for three scenarios. The blue curve represents the brute force approach, with no use of OOVD data. In other words, this is simply a plot of the function . The green curve represents the number of buckets considered when employing OOVD data in a naive way, i.e., with 35 buckets per face. Finally, the orange curve represents the number of buckets in which resulted from processing the refined OOVD data as discussed above.

It can be seen that employing OOVDs results in a significant reduction in the number of buckets that need to be condidered. In fact, for terminals this reduction is of the order over the brute force approach. It also shows that the improved OOVD approach provides at least a factor 12 improvement over the simple OOVD scenario. This is as expected since the reduction from 35 to 3 buckets per face provides nearly a factor 12 in itself.

We now mention a number of statistical properties of optimal -Steiner trees that we observed for trees generated in these experiments.

Firstly, in Figure 6 we have plotted the percentage reduction in network length that is achieved by employing a single Steiner point versus the MST. As is expected, this percentage decreases as increases. In Figure 7 we present a local measure of improvement in length by comparing the length of the Steiner edges added to the MST and the MST edges deleted. Figure 8. The proportion of solutions in which the added and deleted edges are from the same triangle.

We also observed that in many cases the optimal -Steiner tree results from deleting a pair of adjacent edges in the MST. In Figure 8 we have plotted the proportion of times this occurred for our instances. As can be seen, our experiments suggest that this occurs on average between and of the time for uniformly generated instances. Figure 9. The solution for a lattice problem instance contains a cross.

Finally, even though degree- Steiner points are possible in an optimal -Steiner trees, our experiments show that they are rare. In fact, not a single instance we constructed from the random examples contained a degree- Steiner point. Figure 9 shows an example of an optimal tree on a contrived instance which does contain a degree- point.

## References

•  P. K. Agarwal, O. Schwarzkopf, and M. Sharir. The overlay of lower envelopes and its applications. Discrete Comput. Geom., 15(1):1–13, 1996.
•  M. Brazil, C. J. Ras, and D. A. Thomas. A geometric characterisation of the quadratic min-power centre. European J. Oper. Res., 233(1):34–42, 2014.
•  Marcus Brazil, Charl J. Ras, Konrad J. Swanepoel, and Doreen A. Thomas. Generalised -Steiner tree problems in normed planes. Algorithmica, 71(1):66–86, 2015.
•  Maw Shang Chang, Nen Fu Huang, and Chuan Yi Tang. An optimal algorithm for constructing oriented Voronoĭ diagrams and geographic neighborhood graphs. Inform. Process. Lett., 35(5):255–260, 1990.
•  Herbert Edelsbrunner, Leonidas J. Guibas, and Micha Sharir. The upper envelope of piecewise linear functions: algorithms and applications. Discrete Comput. Geom., 4(4):311–336, 1989.
•  Herbert Edelsbrunner and Raimund Seidel. Voronoĭ diagrams and arrangements. Discrete Comput. Geom., 1(1):25–44, 1986.
•  Efi Fogel, Dan Halperin, and Ron Wein. CGAL arrangements and their applications, volume 7 of Geometry and Computing. Springer, Heidelberg, 2012. A step-by-step guide.
•  George Georgakopoulos and Christos H. Papadimitriou. The -Steiner tree problem. J. Algorithms, 8(1):122–130, 1987.
•  Daniel Juhl, David M. Warme, Pawel Winter, and Martin Zachariasen. The GeoSteiner software package for computing Steiner trees in the plane: an updated computational study. Math. Program. Comput., 10(4):487–532, 2018.
•  Clyde Monma and Subhash Suri. Transitions in geometric minimum spanning trees. volume 8, pages 265–293. 1992. ACM Symposium on Computational Geometry (North Conway, NH, 1991).
•  C. J. Ras. Survivable minimum bottleneck networks. Comput. Geom., 50:17–23, 2015.
•  J. H. Rubinstein, D. A. Thomas, and J. F. Weng. Degree-five Steiner points cannot reduce network costs for planar sets. Networks, 22(6):531–537, 1992.
•  Godfried T. Toussaint. The relative neighbourhood graph of a finite planar set. Pattern Recognition, 12(4):261–268, 1980.
•  D. M. Warme, P. Winter, and M. Zachariasen. Exact algorithms for plane Steiner tree problems: a computational study. In Advances in Steiner trees, volume 6 of Comb. Optim., pages 81–116. Kluwer Acad. Publ., Dordrecht, 2000.