Efficient Speed-Up of the Smallest Enclosing Circle Algorithm

1Introduction

The smallest enclosing circle problem is defined as follows. Given a set of 2D 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 brutal force algorithm has the time complexity O(N4). The algorithm tests all possible circles, i.e. all combinations of two points O(N2) and all combinations of three points O(N4). To test if one circle is the smallest enclosing circle takes O(N). This whole brutal force algorithm runs in O(N4) time.

There are two best known algorithms, i.e. the Welzl’s algorithm, which has expected Oexpected(N) running time (Welzl, 1991), and the Megiddo’s algorithm, which has expected Oexpected(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(NlogN) 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+slogh), 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(slogh)≪O(N), thus the expected time complexity of the convex hull algorithm is Oexpected(N) (Skala et al., 2016b).

2Proposed 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 Oexpected(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.

Algorithm 1

Smallest enclosing circle.

2.1Convex 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.

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.

3Experimental Results

In this chapter, we summarize the obtained results of the proposed algorithm for the smallest enclosing circle of points in 2D. 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.

Fig. 2

The smallest enclosing circle for real data sets. The 2D input points are created by projection of all 3D points into xy plane. The models are from the Stanford 3D 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 2D. 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.

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 103 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 102−103 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 {C1,C2,C3,C4}. 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.

Fig. 4

Definition of symbol used in graphs with running times.

Fig. 5

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

Fig. 6

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

Fig. 7

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

Fig. 8

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

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 {C1,C2,C3,C4} 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 {C1,C2,C3,C4}). 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.

Fig. 10

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

4Conclusion

We proposed a simple and efficient speed-up of the Weltzl’s algorithm for computation of the smallest enclosing circle of 2D 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 2D points can be adapted for the computation of the smallest enclosing sphere of 3D points using convex hull computation in 3D (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.