Efficient Speed-Up of the Smallest Enclosing Circle Algorithm

1 Introduction

The smallest enclosing circle problem is defined as follows. Given a set of 2 D points, find a circle with the smallest radius such that all the given points are contained in either inside of this circle or its circumference. The smallest enclosing circle meets all the following requirements:

•

The maximal distance of all points to the center is minimal for the center of the smallest enclosing circle.

•

Given any three points, we can uniquely define a circle, with at least two points of these points on the circumference.

•

The smallest enclosing circle is unique for any set of 2 D points.

•

Given N ⩾ 2 points in the plane, the smallest circle that encloses them contains at least two of the points within its circumference.

The brutal force algorithm has the time complexity O ( N 4 ) . The algorithm tests all possible circles, i.e. all combinations of two points O ( N 2 ) and all combinations of three points O ( N 4 ) . To test if one circle is the smallest enclosing circle takes O ( N ) . This whole brutal force algorithm runs in O ( N 4 ) time.

There are two best known algorithms, i.e. the Welzl’s algorithm, which has expected O expected ( N ) running time (Welzl, 1991), and the Megiddo’s algorithm, which has expected O expected ( N ) running time (Megiddo, 1983). The technique prune and search is used in Megiddo’s algorithm (Megiddo, 1983) to reduce the problem size by removing approximately N / 16 unnecessary points. The algorithm is rather complicated which resulted to a very high multiplicative constant. The reduction needs to solve twice the similar problem where the centre of the enclosing circle is located using the construction of bisectors of points. The randomized incremental construction is used in the Welzl’s algorithm (Welzl, 1991). The algorithm iterates over all points. In case, the point is not inside the actual smallest enclosing circle, the point must be part of the new smallest enclosing circle in the next step. Using this basic knowledge the final minimal enclosing circle is recursively constructed.

The Skyum algorithm (Skyum, 1991) presents a simple iterative O ( N log N ) algorithm for computing the smallest enclosing circle. This is not optimal but its simplicity makes it a better alternative for medium-sized problems than both previous algorithms. The Gau algorithm (Gao and Wang, 2018) computes the minimal enclosing disk in an iterative manner and needs to define some distance parameter delta value that is used in the algorithm.

Other alternative algorithms (Har-Peled and Mazumdar, 2003, 2005) compute the smallest enclosing circle not for all input points but for at least k points. It can also compute an approximation of this circle, which can be useful in some applications as well. The algorithm of the smallest enclosing circle of at least k points is described in Efrat et al. (1993, 1994), too. The paper (Xu et al., 2003) summarizes four different approaches for the location of a minimal enclosing circle for a set of fixed circles.

One of the first sophisticated geometry algorithms developed with many variations of it is the convex hull algorithm. The convex hull is the minimal set of points, that form a boundary for all other points, i.e. all other points are located inside. There are numerous applications for convex hulls: collision avoidance, maximum distance using convex hull diameter for large data sets (Skala, 2013; Skala and Majdisova, 2015), hidden object determination, shape analysis, point inside polygon (Skala and Smolik, 2015).

In this paper, an adapted convex hull algorithm from Skala et al. (2016b), based on Skala (2013), is used to speed-up the computation of the smallest enclosing circle. This convex hull algorithm uses the space subdivision to achieve time complexity of O ( N + s log h ) , where s is number of suspicious points to be on convex hull and h is number of points on convex hull. It is expected that O ( s log h ) ≪ O ( N ) , thus the expected time complexity of the convex hull algorithm is O expected ( N ) (Skala et al., 2016b).

2 Proposed Approach

The smallest enclosing circle contains within its circumference only points from the convex hull. To speed-up the processing time, two step algorithm was developed. The first step is a removal of all points that cannot form the smallest enclosing circle. While the second one is a computation of the smallest enclosing circle using the Weltzl’s algorithm (Welzl, 1991) as it is fast and easy to implement.

In order to obtain a significant speed-up of the input points reduction, the convex hull construction needs to have a lower computational cost than the Weltzl’s algorithm. The convex hull algorithm (Skala et al., 2016b) uses space subdivision technique to speed-up the computation with O expected ( N ) time complexity.

The proposed algorithm for the smallest enclosing circle is summarized in Algorithm 1 and will be described in more detail in the following sub-sections.

OPEN IN VIEWER

Algorithm 1

Smallest enclosing circle.

2.1 Convex Hull Construction

The convex hull algorithm (Skala et al., 2016b) significantly speeds-up the construction of the convex hull by reducing the input points to only a few ones that are suspicious of forming the convex hull. The detailed algorithm is described in Skala et al. (2016b). We summarize only the important steps without details.

OPEN IN VIEWER

Fig. 1

Convex hull construction. Blue lines represent the convex hull.

The first step is to find an estimation of the axis-aligned bounding box (AABB) if not known. This is done by searching only 10 % of input points. A convex polygon is estimated from those points defining the estimated AABB. All points that are inside of this polygon cannot form the convex hull and thus are discarded (the gray part in Fig. 1).

The next step is to create a star shape division of the remaining points. In each cell, we determine points with the maximal distance from the centre, that can form the convex hull (there can be more points in each cell). All other points are discarded. The remaining points (connected with red lines in Fig. 1) are the result of the reduction. Now, the actual algorithm for convex hull construction using only the suspicious points is used.

2.2 Smallest Enclosing Circle Computation

The Welzl’s algorithm for the smallest enclosing circle is recursive. Originally, the algorithm is randomized. However, in this approach, it is used without randomization as the input points are already sorted in the required order enabling speed-up of the computation significantly.

One limitation of the Weltzl’s algorithm is the recursion as it allows to process only limited number of input points, due to the depth of recursion. We overcome this limitation by reducing the input points to only points on the convex hull. Now, the amount of points to be processed is limited only by available RAM memory.

The original algorithm selects points totally in a random manner. However, if the first selected points form a circle, which is big enough to contain almost all the points, then the algorithm will speed-up. In the first step, we locate the two points { C 1 , C 2 } with maximum distance from each other using the algorithm (Skala and Majdisova, 2015). Next, we locate the point C 3 , which has the maximal distance from the center of previously located two points { C 1 , C 2 } . Last, we locate one more point C 4 , which has the maximal distance from point C 3 . All of those four points are moved to the end of input points for modified Weltzl’s algorithm and the rest of the points are randomly shuffled. In the case, when the smallest enclosing circle is defined only by two points, then these two points are exactly { C 1 , C 2 } and the points { C 3 , C 4 } are the same as { C 1 , C 2 } .

The input for the Weltzl’s algorithm is a set of points P. In each step, the algorithm selects the last one point p from P, and recursively finds the smallest circle containing P − { p } , i.e. all of the other points in P except of point p. If the point p is also included in the returned circle, then it is the smallest enclosing circle for the whole set of points P.

Otherwise, if the point p lies outside the circle, it must lie on the boundary of the resulting circle. In the next step, it recurses with p as an additional point in R (points known to be on the boundary for already tested points from P).

The recursion terminates when P is empty, and a solution can be found from the points in R. In case of 0 or 1 points, the circle is only none or one point. The smallest enclosing circle for 2 points is defined, that has its centre in the middle of the two points and radius half of the distance between the two points. The last case are 3 points, where the smallest enclosing circle is the circumcircle of the triangle described by the points.

When R contains 3 points, the recursion terminates, as all points in P are already inside of the circle formed by R.

3 Experimental Results

In this chapter, we summarize the obtained results of the proposed algorithm for the smallest enclosing circle of points in 2 D . The proposed approach is capable to compute the smallest enclosing circle for both synthetic and real data sets as well. The synthetic data sets were used to measure the performance of the proposed algorithm. The result of the proposed approach on real data sets is visualized in Fig. 2. The 3 points that define the smallest enclosing circle for bunny were actually selected between the first 4 points that are processed at the beginning of the Weltzl’s algorithm. For the dragon data set were 2 out of 3 points that define the smallest enclosing circle selected to be in the first 4 points that are processed by Weltzl’s algorithm.

OPEN IN VIEWER

Fig. 2

The smallest enclosing circle for real data sets. The 2 D input points are created by projection of all 3 D points into x y plane. The models are from the Stanford 3 D scanning repository (http://www.graphics.stanford.edu/data/3Dscanrep/).

To test the proposed approach, we used several data sets with different types of point distribution in 2 D . Some of the distributions are well known, randomly distributed uniform points inside a unit square or inside a unit circle and points with Gaussian distribution. Then we also used Halton points and Gauss Ring points. The last two distributions are described in Skala et al. (2016b). The visualization of all tested distributions of points together with their convex hulls is presented in Fig. 3.

OPEN IN VIEWER

Fig. 3

Distributions of points used for testing. The blue line represents the convex hull of the data set, i.e. the initial points for Weltzl’s algorithm.

The main purpose of the proposed approach is the speed-up of the Weltzl’s algorithm. The running time of the algorithm for the smallest enclosing circle depends on the distribution of points and of course on the number of input points. We performed 10 3 tests for each points distribution from Fig. 3 and for each different number of input points. The running time of the algorithm for one input data set was measured as running time of 10 2 − 10 3 repetitions divided by the number of repetitions. The random generator was initialized with the same seed number for each of those repetitions. The times were measured for the original Weltzl’s algorithm and in the first phase also for the proposed approach without the selection of four best candidates { C 1 , C 2 , C 3 , C 4 } . As the algorithm is randomized, the running time depends on the randomization and thus we computed the minimal and maximal running time for each configuration as well as the 33th% fastest and 33th% slowest time. The visualization of these four running times is described in Fig. 4. The resulting running times are visualized in Figs. 5–9. It should be noted, that there is a logarithmic scaling on both axes. The horizontal axis on graphs with proposed approach presents higher maximal number of input points compared to graphs with standard Weltzl’s algorithm.

OPEN IN VIEWER

Fig. 4

Definition of symbol used in graphs with running times.

OPEN IN VIEWER

Fig. 5

Running times in [ μ s ] (vertical axis) for different number of uniform points in circle (horizontal axis).

OPEN IN VIEWER

Fig. 6

Running times in [ μ s ] (vertical axis) for different number of Gauss Ring points (horizontal axis).

OPEN IN VIEWER

Fig. 7

Running times in [ μ s ] (vertical axis) for different number of Gauss points (horizontal axis).

OPEN IN VIEWER

Fig. 8

Running times in [ μ s ] (vertical axis) for different number of Halton points (horizontal axis).

OPEN IN VIEWER

Fig. 9

Running times in [ μ s ] (vertical axis) for different number of uniform points in square (horizontal axis).

It can be seen that for the original Weltzl’s algorithm there is a difference between the fastest and the slowest running time for one distribution and one number of points almost every time at least 10 times different. This is a huge difference in the algorithm behaviour. The proposed approach without the selection of four best candidates { C 1 , C 2 , C 3 , C 4 } has the difference between slowest and fastest time mostly around 3.2 times, i.e. much better ratio resulting to stability of algorithm behaviour.

It can be also seen that we measured the running times of Weltzl’s algorithm for a lower maximal number of input points. The reason for this is that this was the maximal size of the input data set, as there is a significant limitation due to the maximal depth of recursion.

In the second set of tests, we measured the running times with the same distributions and the same number of points for our proposed approach (with all steps of the algorithm, i.e. also with the selection of four best candidates { C 1 , C 2 , C 3 , C 4 } ). The fastest times are the same as in Figs. 5b–9b but the slowest times are lower. The slowest times are now only about 1.9 times slower than the fastest times. This is a great improvement from the previous 3.2 times difference.

We also computed the speed-up of our algorithm to the Weltzl’s algorithm. We used the average running times to compute the speed-up (see Fig. 10). It can be seen that with an increasing number of input points, the speed-up increases as well. This proves, that our proposed approach is not only faster, but it has also better time complexity than the original Weltzl’s algorithm.

OPEN IN VIEWER

Fig. 10

Speed-up of average running times of the proposed method compared to the standard Weltzl’s algorithm.

We proposed a simple and efficient speed-up of the Weltzl’s algorithm for computation of the smallest enclosing circle of 2 D points. The proposed approach is easy to implement. It reduces the randomness of the algorithm and speeds-up the computation, too. The proposed algorithm also improves the computational complexity, which is beneficial for a higher number of processed input points. Also, the maximal number of points, that can be processed, is significantly increased, as the required depth of recursion is dramatically decreased.

The proposed algorithm for computation of the smallest enclosing circle of 2 D points can be adapted for the computation of the smallest enclosing sphere of 3 D points using convex hull computation in 3 D (Skala et al., 2016a). Another possible improvement could be usage of the smallest enclosing circle/sphere of the points forming axis-aligned bounding box as the initial step because all points inside cannot form the smallest enclosing circle/sphere.