The Earth is nearly flat: Precise and approximate algorithms for detecting vulnerable regions of networks in the plane and on the sphere

Several recent works shed light on the vulnerability of networks against regional failures, which are failures of multiple pieces of equipment in a geographical region as a result of a natural disaster. To enhance the preparedness of a given network to natural disasters, regional failures and associated Shared Risk Link Groups (SRLGs) should be first identified. For simplicity, most of the previous works assume the network is embedded on a Euclidean plane. Nevertheless, they are on the Earth's surface; this assumption causes distortion. In this work, we generalize some of the related results on the plane to the sphere. In particular, we focus on algorithms for listing SRLGs as a result of regional failures of circular or other fixed shape.


INTRODUCTION
Serious network outages are happening with increasing frequency due to disasters (such as earthquakes, hurricanes, tsunamis, tornadoes, etc.) that take down almost every equipment in a geographical area (see [9] for a recent survey conducted within COST Action RECODIS [21] on strategies to protect networks against large-scale natural disasters). Such failures are called regional failures and can have many locations, shapes, and sizes.
Due to the huge importance of telecommunication services, improving the preparedness of networks to regional failures is becoming a key issue [6, 8, 10-12, 16, 17, 19, 24]. Roughly speaking, protecting networks against regional failures is dealt with either by using geometric tools [1,7,19,[25][26][27]29] or by reducing massively the problem space to a set candidate locations of failures [6,10,11,16,17,24]. Nevertheless, both approaches require a detailed knowledge of the geometry of the network topology, such as the precise GPS coordinates of nodes and cable conduits' routes, and the statistics of past disasters.
In many works, regional failures are computed by transforming the geographical coordinates of an existing network into a plane, which introduces distortion. Depending both on the geographical area of the network and on the transforming procedure, network is left intact. Among the possible geometric failure shapes, the most natural one is the circular disk, as it is compact, and is invariant to rotation. One possibility is to overestimate the possible failures with circular disks (or any other fixed geometric shape), which yields short SRLG lists. However, it is not clear what is the cost of this overestimation. In most of this work, we choose to overestimate the disasters by circular disks with a maximum size according to one among many possible measures, while we also briefly tackle overestimations with different geometric shapes.
When talking about (maximal regional) disk failures, the most natural measure is the disk radius, which represents the maximum geographical coverage of the natural disaster. Nevertheless, since the network density is usually not homogeneous (i.e., there are more nodes and links in crowded geographical areas than in non-crowded areas) the number of network elements (either nodes or links) contained by the disk are also two useful measures (it is natural that more SRLGs are needed in crowded areas and less crowded areas can be covered with fewer SRLGs). Therefore, in most of this work, we will concentrate on the following three types of SRLG lists: • maximal r-range SRLG list: list of maximal link sets that can be hit by a disk with radius at most r.
• maximal k-node SRLG list: list of maximal link sets that can be hit by a disk hitting at most k nodes.
• maximal k-link SRLG list: list of maximal link sets that can be hit by a disk hitting at most k links.
To distinguish between lists obtained from planar and spherical representation, we will include attribute planar or spherical in the list names (e.g., maximal spherical r-range SRLG list) when needed for clarification.
It turns out that in all three mentioned cases, the size of the maximal SRLG list is linear in the network size in practice, and can be computed in low-polynomial time both in the planar case (maximal r-range list: [25], k-node: [29], k-link: [20]) and the spherical case.
This paper is an extension of [30], which is, to the best of our knowledge, the first study on enumerating regional SRLG lists in a spherical model. Also, in [30], the first comparison of the spherical and planar representations was provided. It is shown through examples that precise polynomial algorithms could be designed for spherical representation of the networks. In our experience, these algorithms are only 2 times slower than their planar pairs. We also believe that using our approach, resilient routing and network design results [14] can be further enhanced. Compared to the conference version of our paper [30], we have (1) included Section 4.2 presenting approximate algorithms on enumerating the maximal failures induced by arbitrary disaster shapes, (2) provided an enhanced mathematical analysis of algorithms presented Sections 4.1 and 4.2, and (3) in Section 5 added simulation results showing that the difference between the planar and spherical representations of the network can result in different SRLG lists even in case of networks having a geographic extension as small as 100km.
As there are many mathematical derivations in the rest of the paper, we would like to summarize the concepts in plain text once again here for the sake of readability. As learned from previous studies, all of r-range, k-node, and l-link lists can be precisely calculated in low-polynomial time of the network size in case of planar representation of the network. Our first goal was to show that considering spherical embeddings of the networks the possibility of designing low-polynomial time algorithms for determining these SRLG lists remains. We demonstrated this phenomenon in Section 3 by designing an algorithm capable of determining the r-range SLRG list both in the planar and spherical case in low-polynomial time of the network elements. The existence of fast precise algorithms is good news, however, their drawback is that intuitively the faster they are, the harder they are to implement. This fact motivates our second goal, namely designing a framework of simple and fast algorithms capable of determining all the mentioned SRLG lists in both planar and spherical representation with enough precision, which are presented in Section 4.1. Our third and final goal is to show how simple approaches can be applied to more general models too: while in most of this study, we concentrated on disasters having a shape overestimated by a disk, Section 4.2 gives an outlook for designing easy-to-implement algorithms for essentially arbitrary disaster shapes.
The remaining of the paper is organized as follows. Section 2 describes the network representation model together with the assumptions made. In Section 3, we present an example of a polynomial algorithm for computing maximal SRLG lists, handling both the planar and spherical cases. While in Section 4.1, a faster and more flexible approximate approach is presented for solving the same problem, Section 4.2 gives an outlook for designing algorithms for arbitrary disaster shapes. Simulation results are presented in Section 5 and, finally, we draw the conclusions in Section 6.

MODEL AND ASSUMPTIONS
Throughout the paper, we will consider two types of embeddings of the network: embedding in Euclidean planar and spherical geometry.
The network is modeled as an undirected connected geometric graph  = (V, ) with n = | V | ≥ 3 nodes and m = | | edges stored in a lexicographically sorted list. The nodes of the graph are embedded as points in the Euclidean plane or sphere, and their precise coordinates are considered to be given in 2D and 3D Cartesian coordinate system in the planar and spherical cases, respectively. Note that if coordinates are given in polar system (in case of spherical geometry), one can easily transform them to Cartesian at the very beginning.
When speaking of planar geometry, for each edge e, there is a polygonal chain (or simply polyline) e l in the plane in which the edge lies (see Figure 2). Parameter will be used to indicate the maximum number of line segments a polyline e l can have. Naturally, in spherical case, the polyline of an edge refers to a series of geodesics. Note that this model covers special cases when edges are considered as line segments (geodesics).
It will be assumed that basic arithmetic functions (+, −, ×, ∕, √ ) have constant computational complexity. For simplicity, we assume that nodes of V and the corner points of the containing polygons defining the possible route of the edges are all situated in general positions of the plane, that is, there are no three such points on the same line and no four points on the same circle, and in the spherical case there are no antipodal nodes or breakpoints and no great circles of geodesics of polylines cross the North pole. In this study, our goal is to generate a set of SRLGs, where each SRLG is a set of edges. Note that from the viewpoint of connectivity, listing failed nodes besides listing failed edges has no additional information. We consider SRLGs that represent worst-case scenarios the network must be prepared for and, thus, there is no SRLG which is a subset of another SRLG.

Model for circular disk-shaped disasters
In most of this study, it will be assumed that disasters are either having a shape of a circular disk or they are overestimated by a circular disk.
We will often refer to circular disks simply as disks. The disk failure model is adopted, which assumes that all network elements that intersect the interior of a circle c are failed, and all other network elements are untouched. We emphasize that in this model when we say e is hit by c, it does not necessarily mean that e is destroyed indeed by c; instead, it means that there is a positive chance for e being in the destroyed area. In other words, this modeling technique does not assume that the failed region has a shape of a disk, but overestimates the size of the failed region in order to have a tractable problem space.
Definition. Let  p and  s denote the set of all disks in the plane and the set of all disks on the sphere, respectively. For both geometry types g ∈ {p, s}, let  g r ,  g k and  g l denote the set of disks part of  g having radius at most r, hitting at most k nodes of V and hitting at most l links of , respectively. Based on the above definition, we define the set of failure states that a network may face after a disk failure, with a maximal measure.
Definition. For all geometry types g ∈ {p, s} and SRLG type t ∈ {r, k, l}, let set F( g t ) denote the set of edges which can be hit by a disk c ∈   Table 1 gives an overview on the corresponding notations and denominations of the SRLG list types we focus on this paper.  One aim of this study is to propose fast algorithms computing these lists for various sizes of m.

Model for disasters with arbitrary shape
In Section 4.2, we give an overview of simple algorithms handling disasters which have a shape of a fixed non-self-intersecting closed polytope in the plane having boundaries consisting of line segments and arcs. The disaster can occur in any physical area and with any orientation. A link is hit is it has at least a common point with the disaster. The model is basically the same as the one described in the previous section, the only difference being the disaster shape.

PRECISE ALGORITHMIC APPROACHES FOR ENUMERATING MAXIMAL FAILURES
As mentioned before, determining the planar lists M p r , M p k , M p l is relatively well studied. It remains a question of how much the distortion of maps can affect the calculated SRLG lists. The answer is that it heavily depends on the projection used to make the map. For example, while the stereographic projection affects significantly the distances, but in contrast to many other projections it has the nice property of mapping spherical disks to planar disks [22] (fact also used in Appendix). One approach for calculating spherical lists would be to adapt existing algorithms to spherical geometry demonstrating the interoperability between these geometries. However, in this paper, we follow an approach simpler to present and avoiding trigonometric calculations via applying the projection in both directions for numerous times. In other words, some steps of the algorithm are performed on the plane, while others on the sphere.
In the following, we extend a precise algorithm for determining M p r (see [25]) to an algorithm computing M p r or M s r depending on the geometry of the input. In the rest of this section, we present this extended algorithm.

Smallest enclosing disks
Let us make the following definition for the sake of clarifying the intuition.
Definition. Let a disk c be smaller than disk c 0 , if c has a smaller radius than c 0 , or if they have an equal radius and the center point of c is lexicographically smaller than the center point of c 0 . Among a set of circles S c , let c be the smallest if it is smaller than any other circle in S c .
Definition. Let F ⊆ E be a finite nonempty set of edges (not necessarily a failure). We denote the smallest disk among the disks enclosing the polylines of F by c F and we say c F is the smallest enclosing disk of F.
It is not difficult to see that c F always exists for line segments or geodesics (depending on the geometry), and thus, by mapping the corresponding segments/geodesics together we can deduce that the definition is correct for polylines too. The key idea of our approach is that we can limit our focus only on the smallest enclosing disks c F . The consequence of the next proposition is that the number of smallest enclosing disks c F is not too large.

Proposition 1. Let H be a nonempty set of polylines of edges with smallest enclosing disk c H . Then there exists a subset
Definition. Let  denote the set of maximal edge sets hit by a smallest enclosing disk.
By Proposition 1 we have:

Lemma 3. Let H be a set of line segments in the plane or geodesics on the sphere, |H | ≤ 3. Then c H can be determined in O(1) time.
The proof of the Lemma is relegated to the Appendix.

Polynomial algorithm for determining maximal failures
In this section, we repeat an extension of the basic algorithm provided by [25] which handles both spherical and planar inputs. There are two key facts inspiring this algorithm. Firstly, based on Proposition 1:

On improved complexity bounds
The results from the previous section can be easily improved using parametrization and some computational geometric tricks.
The first observation is that for any meaningful radius of the disk failure most of the network will remain intact. However, failures with the same radius taking place in a crowded area tend to take down more equipment than the ones in sparsely inhabited areas. This motivates the introduction of a graph density parameter: Definition. For every r ∈ R + 0 , let r be the maximum number of edges which can be destroyed by a disk with radius at most r.
This r is considered to be small in case of small r values. Another observation is that there are not many more network edges than nodes. This is formalized in the upcoming Claim 1.
Informally speaking, we denote the set of crossing points of the edges by X. A more formal definition follows.
Definition. Let X be the set of points P in the plane on which no node element of V lies and there exist at least 2 edges which have polylines having a finite number of common points crossing each other in P. Let x = | X|.
Despite the fact that on arbitrary graphs x can be even O(n 4 ), in backbone network topologies typically x ≪ n because a node is usually installed if two cables are crossing each other. This gives us the intuition that G is "almost" planar, and thus it has few edges. Claim 1. The number of edges in G is O(n + x). More precisely for n ≥ 3 we have m ≤ 3n + x − 6.
Proof. If  is embedded in the plane, do the following. Let G 0 (V ∪ X, E 0 ) be the planar graph obtained from dividing the polylines of edges of G at the crossings. Since every crossing enlarges the number of edges at least by two, |E 0 | ≥ m + 2x. On the other hand, If  is embedded in the sphere, we can project it to the plane with stereographic projection, repeat the former arguments then apply an inverse projection to the sphere. time.
Corollary 9 proposes that M g r can be determined in quartic time of n in practice. On the other hand, Algorithm 1 has its limits of speed: because of the three nested for-loops, it runs in Ω(n 3 ) time. In order to achieve better results, the algorithm would have to be changed. For the planar case, Reference [25] gives an algorithm which runs in O((n + x) 2 5 r ) time for = 1 (i.e., the edges are considered as line segments there). Furthermore, we are convinced that an algorithm with parametrized running time near-linear in network size could be achieved for determining M g r (and also for determining M g k and M g l , despite the fact that they can be computed based on very different theories). However, presenting this kind of algorithm would exceed the limits of this paper.

APPROXIMATE ALGORITHMS AND IMPLEMENTATION ISSUES
It is always good to have fast precise polynomial algorithms for solving a given problem. However, this approach also has some disadvantages: (1) the lower complexity a precise algorithm for determining a maximal circular SRLG list has, the harder to implement and prove its correctness and complexity; (2) designing algorithms for computing different types of maximal SRLG lists need totally different mathematics. Moreover, in most cases, the available geographical data of networks are inaccurate. Adding this fact to the inconveniences of the precise approach results into the idea of designing some approximate algorithms that are able to compute these lists with enough precision.
The main idea behind this class of algorithms that instead of keeping the original shape and size of the disaster area and taking in count all the infinitely many possible disaster centers, we slightly overestimate the disaster, letting us detect all the same (or occasionally a bit larger) hit link sets while going through only a finite number of centers.

Approximate algorithms for circular disk failures
In this section, we present an approximate approach suitable for computing all types of maximal SRLG lists defined in Section 2.
Definition. For a point P (in the plane or on the sphere) and node v ∈ V, let the node-distance couple be [v, d(v, P)], where d(v, P) is the distance between v and P. Let e(P) be the list consisting of the link-distance pairs of all links e ∈ , sorted according to the lexicographical order of the links. Let e(P) hit be the sorted list of links not further from P than r.

Proposition 10. For a given point P, both e(P) and e(P) hit can be computed in O((n + x) ) time.
Clearly, both node-distance lists and edge-distance lists can be determined quickly. The plan is to determine these lists for enough points which are also placed well enough to be able to determine the maximal SRLG lists based on these node-distance and edge-distance lists.
Definition. Let  denote the set of points P for which we want to construct the link-distance lists e(P).
Let us stick to planar geometry for a moment. Intuitively, we can calculate M p r by including the grid points of a sufficiently fine grid (let us say containing 1 km × 1 km squares) in . On the sphere, we should choose a similar nice covering. It is possible that we have some extra short links, thus for calculating the k-node and k-link list we should include some extra points in . For example, by adding some random points of polylines of each edge and some point near every node, we can solve this issue.
Definition. Let  be the maximal distance of any geometric location from the (closed) convex hull of the geometric embedding of graph G to the closest point of set , that is,  ≔ max t∈conv(G) min p∈ dist(p, t).
Definition. Taking two sets E 1 and E 2 , we denote the relationship of the sets with E 1 ⊒ E 2 if and only if for all e 2 ∈ E 2 there exists an e 1 ∈ E 1 , such that e 1 ⊇ e 2 . Based on Theorem 11, if one wants to protect disasters caused by disks with radius r, it is only necessary to run Algorithm 3 initializing the radius as r +  .
Comparing Corollaries 13 and 9, we can see that despite the fact that the approximate Algorithm 3 is much simpler to implement, and taking into account that disasters are not precisely circular and the chosen radius is arbitrary, it clearly outperforms the precise Algorithm 1.

Approximate algorithms for disasters with arbitrary shape
Understanding how to deal with circular disk failures is a good start; however, one should consider other disaster shapes too. In [28], it is proven that, in a similar problem setting (the model of [28] considers that a link is hit by a disaster exactly if at least one of its endpoints is inside the disaster region. We believe our failure model is more realistic), there are a polynomial number of maximal failures caused by disasters having elliptic or polygonal (e.g., rectangular or equilateral triangular) shape. Again, as engineering fast precise algorithms for determining SRLG lists similar to M r , M k or M l but for arbitrary disaster shapes instead of a circle is not trivial, we are going to discover the possibilities of approximate algorithms similar to the one described for determining M r . In short, while the disk is invariant to rotation, now one should consider also the different orientations of the fixed shape. We shall discretize the continuous rotation of the shapes via taking only a possible orientations of the shape: Definition. Let F r be a set of non-self-intersecting closed polytopes embedded in the plane having boundaries consisting of line segments and arcs, for which (1) every f 1 , f 2 ∈ F r , f 1 can be moved into f 2 via a translation and a rotation, and (2) the smallest fitting disk of any f ∈ F r has a radius of r.
Let T be an arbitrary point of a polytope f ∈ F r . Let N(F r ) = N  + a,T (F r ) be the consisting of the elements f ′ of N  (F r ), each f ′ extended as f ′ ′ with the smallest d a, T -neighborhood of f ′ (where distance d a, T is a function of the extended shape f ′ , the number of orientations of the shape considered a, and the center of rotation T ∈ intf ) in such a way that: ⋃ By the former definition, we have the following. For any possible place and orientation of a disaster with shape f , in which f intersects at least a network element, there exists a p ∈  and an orientation of the extended disaster shape f ′ ′ for which the area destroyed by f ′ ′ contains the area destroyed by f . This also holds if the orientation is taken from among the a investigated orientations. Thus the edge set hit by f ′ ′ is a superset of the edge set hit by f .
Definition. Let denote the number of line segments and arcs needed to describe a shape from F r . Let r be the maximal number of edges an f ∈ F r can hit. Proof. Clearly, the algorithm computes M g N(F r ) . Based on Proposition 14 proposing that as  → 0 and a → ∞, the necessary inflation of the original shape f tends to 0, we also can see that for any given network, lim  → 0 Regarding its complexity, we can argue as follows.
To calculate  , we observe that all the points in N r (G) (the smallest r-neighborhood of the convex hull of G) maximizing the distance from  lie on the intersection of at least 3 Voronoi cells computed on point set , restricted to N r (G), and the outer region of N r (G) taken as an extra cell. The Voronoi regions of  can be computed in O(|| log ||) time both in the plane and on the sphere (plane: Section 7.2 in [5], sphere: special case of Corollary 2 in [18]). Using Claim 1, one can see that the number of line segments and arcs (or geodesics on the sphere) needed to describe N r (G) is O ((n + x) ), thus we claim the total complexity of computation and description the modified Voronoi diagram described before is O((n + x) || log ||). Based on the graph representation of the diagram one can determine  in the proposed complexity.
We claim that Line 2 can be done in O( ) time if f is convex. In case of non-convex failure shapes, f calculations can become more complicated (e.g., holes in f can disappear in f ′ ′ ), but clearly, f ′ ′ can be determined in time polynomial in in this case too.
For a given P ∈ , translation of f (Line 4) can be done in O( ) time.
For a given P ∈  and z ∈ { 0, 2 a , … , 2(a−1) a } , rotation (Line 6) can be done in O( ) time. For a given point and direction, Line 7 needs O((n + x) + r ) operations, as follows. For every edge e ∈  one has to check whether its polyline e l intersects the boundary of f z or if not, is e l situated entirely in f z , both can be checked in O( ). Then, refreshing M g N(F r ) with the edge set hit by f z can be done in O( r ) time, since edges are stored in ordered lists.
We can conclude that Algorithm 4 is correct and runs in the proposed complexity. ▪ Algorithm 4 clearly runs in time polynomial in the input size (even for non-convex disaster shapes f ), and can be applied for a wide variety of disaster shapes.
Corollary 17 will offer a more intuitive complexity result on Algorithm 4: ) .
If we compare Corollaries 17 and 13, we can see that unless the shape f of disasters is very complicated, the only real additional complexity compared to the circular disk failure case arises from the fact that that usually f is not invariant to rotation.
Comparing Corollaries 17 and 9, we can see that although the approximate Algorithm 3 handles the much more complex problem induced by disasters having an arbitrary given shape f , it has a lower complexity compared to the precise Algorithm 1 dealing with only circular disasters, showing the strength of the proposed framework of approximate algorithms.

SIMULATION RESULTS
In this section, we present numerical results that validate our approximate approach for circular disk failures presented in Section 4.1, and demonstrate the use of the proposed algorithms on some realistic physical networks. The algorithm was implemented in Python 3.5 using various libraries. Distance functions were implemented from scratch. No special efforts were made to make the algorithm space or time optimal. Run-times were measured on a commodity laptop with Core i5 CPU at 2.3 GHz with 8 GB of RAM. The output of the algorithm is a list of SRLGs such that no SRLG contains the other.
We interpret the input topologies in two ways: polygon, where links are polygonal chains, and line, where the corner points of the polygonal links are substituted with nodes (of degree 2). Here links are line segments.

Extreme geographical extension makes difference
We found that running times for spherical representations were ∼2 times slower than the planar ones in case of most networks (see Table 2). The only exception is when the network has an extreme geographic extension (e.g., AboveNet), in this case, the obtained SRLG lists tend to be longer (Figure 3 demonstrates this in case of k-link lists) causing a slight increase both in the parameter and in the running time.
Another issue which can be noticed is related to the achievable preciseness using the approximate approach. Based on Theorem 11, running time is proportional to ||; given this and the running times collected in Table 2, we can deduce that if the drop of price of computation power remains for an additional short time period, one will be able to run these simulations As a rule of thumb, we can deduce that the difference between the planar and spherical representations of the network can result in different SRLG lists even in the case of networks having a geographic extension as small as 100km.

CONCLUSION
We investigated the problem of generating SRLG lists of networks. We found that the known precise low-polynomial SRLG generating techniques can be modified in order to fit the spherical geometry, allowing us to generate SRLG lists with more precision. A framework of easy-to-implement approximating algorithms for determining the SRLG lists in both planar and spherical representation was also presented.
In our experience, SRLG lists generated using the spherical representation of the networks are different from the planar ones, and also they tend to be longer, especially in case of extreme geographical extension. In the case of our implementation, enumerating SRLG lists in case of spherical representations was typically 2 times slower than in the planar case. The difference between the planar and spherical representations of the network can result in different SRLG lists even in the case of networks having a geographic extension as small as 100km.
While in most of our study, we supposed the disasters destroy the network in an area of a circular disk, some results were generalized to essentially arbitrary disaster shapes. For planar geometry, this problem is already solved, see Theorem 3 of paper [25]. It remains to prove it in case of spherical embedding.
Let e 1 , e 2 , e 3 be three geodesics on the sphere. The endpoints (p i1 , p i2 ) are given by Cartesian coordinates (x i1 , y i1 , z i1 ), (x i2 , y i2 , z i2 ). Let p i3 be an arbitrary point inside e i . Points p i1 , p i2 and p i3 determine the great circle on the sphere containing geodesic e i .
We will project geometric objects on the sphere to the plane using the stereographic projection from the north pole, which has the property that the image of a spherical circle will be a circle on the plane, or in special case, if it contains the north pole, its image is a line [22]. Note that for the sake of simplicity it is assumed that no great circle investigated crosses the north pole.
Projecting the 9 spherical points onto the plane, we receive q ij points given by Cartesian coordinates (x ij , y ij , z ij ) → . We denote the images of e 1 , e 2 , e 3 by arcs f 1 , f 2 , f 3 . Calculating the radius and center point of the containing circle c i for arc f i requires constant number of coordinate geometric steps. Let (x 1 , y 1 ), (x 2 , y 2 ), (x 3 , y 3 ) and r 1 , r 2 , r 3 be the Cartesian coordinates of the center of containing circles, and radii, respectively. The smallest enclosing disk c H on the sphere has an image c ′ H on the plane. However, the parameters of c ′ H are different from the parameters of c H , and the images of the fitting points of c H and e i are the fitting points of c ′ H and f i . That inspires the plan to find all the fitting circles of f i (i.e., those which have exactly one common point with each f i or which have one common point with two of them and containing the third) on the plane, project them back onto sphere and select the smallest among them, as that is the smallest enclosing disk of e 1 , e 2 , e 3 . Thus, we need to find the potential best fitting circles in the plane.
It is possible that the disk fits for two arcs and includes some points of the third. We can choose two arbitrary arcs in 3 ways. Choosing f 1 , f 2 we must calculate the distance of the two arcs and use it as the diameter of the potential disk. On each arc, the distance is determined by an inside point or one of the boundary points. Calculating the distance of two points, a point and a circle or two circles have both constant complexity. So, in this case, 3 ⋅ 3 2 ⋅ O(1) calculations are required.
If the smallest disk touches all of the arcs there are also more different cases. Each arc can be touched on a boundary point or on an inside point (3 3 cases). Fortunately, fitting a circle is already solved in all of the cases and these are known as the problems of Apollonius [4].
If the smallest enclosing disk touches all three arcs f 1 , f 2 , and f 3 , we have three cases for each arc f i : the disk either touches the arc in an interior point or at one of its endpoints. In the former case let (x 1 , y 1 ) and r i be the Cartesian coordinates of the center point and radius of the containing circle of arc f i , respectively. In the latter case, let (x 1 , y 1 ) be coordinates of the endpoint itself, while let r i be 0. Numbers s 1 , s 2 , and s 3 are ±1 representing that the fitting circle touches on the outside or on the inside of the containing circles of c 1 , c 2 , and c 3 (2 3 different possibilities to be checked in each case). Parameters x s , y s , and r s of the fitting C circle can be calculated by solving the following system of equations [2]: (x s − x 1 ) 2 + (y s − y 1 ) 2 = (r s − s 1 ⋅ r 1 ) 2 (x s − x 2 ) 2 + (y s − y 2 ) 2 = (r s − s 2 ⋅ r 2 ) 2 (x s − x 3 ) 2 + (y s − y 3 ) 2 = (r s − s 3 ⋅ r 3 ) 2 .
The system is quadratic, thus it can be solved in a constant number of arithmetic calculations. The complexity of these calculations all together are 3 3 ⋅ 2 3 ⋅ O(1).
Finding the minimal radius of potential disks requires 3 5 − 1 comparisons between diameters. Using this method the minimal disk for e 1 , e 2 , e 3 can be determined in O(1) time. However, the algorithm could be improved by using pre-computations for the edges, excluding some possible disks already on the plane instead of transforming back or fixing s i in case of boundary points. Note that only a constant number of basic arithmetic functions (+, −, ×, ∕, √ ) were used during the computation. ▪