Approximate minimal bounding sphere

Problem statement

Suppose you have a large number of points in 3D space and you want to find a sphere containing all the points. You’d like the sphere to be as small as possible, but you’re willing to accept a slightly larger sphere in exchange for simplicity or efficiency.

False starts

If you knew the center of the sphere, the solution would be easy: find the distance from every point to the center, and the large distance is the radius. But you don’t know the center of the bounding sphere when you start.

What if you found the two points the furthest apart from each other? Maybe that would give you the diameter of the bounding sphere, and the center would be the midpoint of the line connecting these two points.

There are a couple problems with this approach. First, how would you find the two points furthest from each other? The most obvious approach would be to compute the distance between all pairs of points. But if you have n points, you have n(n+1)/2 pairs of points. That would make the algorithm O(n²).

If you do know the two points that are furthest apart, that gives you a lower bound on the diameter of the bounding sphere, but not the diameter itself. For example, suppose we’re working in two dimensions and discover that the longest line between two points connects (-10, 0) and (10, 0). You might try centering a circle of radius 10 at the origin. But what if there’s a third point (0, 11)? It doesn’t change which two points are maximally distance, but it’s a distance of 11 from the origin.

Jack Ritter’s algorithm

A fairly simple algorithm by Jack Ritter runs in linear time, i.e. O(n), and is guaranteed to produce a sphere containing the given points. But the sphere it returns might be a little bigger than necessary. Ritter says [1] “the sphere calculated is about 5% bigger than the ideal minimum-radius sphere.”

When I read that, my first thought was “What kind of distribution is he assuming on points to come up with that 5% estimate?” The expected efficiency of the algorithm must depend on the the distribution on the input points.

Ritter gives two examples. The first started with a sphere “randomly populated with 100,000 points.” Random according to what distribution? I assume uniform. The calculated sphere had a radius 4% larger than optimal. The second example was a cube “randomly populated with 10,000 points.” Again I assume the distribution was uniform. The calculated sphere had a radius that was 7% larger than optimal.

I’ve seen a couple sources that say Ritter’s algorithm is between 5% and 20% from optimal. I’m curious what the 20% figure is based on; the examples below fit much better than that. Maybe the fit is worse for point sets with less symmetry. The number of points may also effect the expected efficiency.

Testing Ritter’s algorithm

To test Ritter’s algorithm, I grabbed an implementation in Python from here. (There’s a flaw in this implementation. See update below.)

Cube example

First, I took 10,000 uniform samples from the unit cube. I started with this example since sampling from a cube is easy.

I changed the last four of the random samples to be corners of the cube so I’d know the exact answer.

    from scipy.stats import uniform
    x = uniform.rvs(size=(10000, 3)).tolist()
    x[-1] = [0, 0, 0]
    x[-2] = [0, 0, 1]
    x[-3] = [1, 0, 0]
    x[-4] = [1, 1, 1]
    print( bounding_sphere(x) )

This returned a sphere of radius 0.8772 centered at (0.5010, 0.4830, 0.4968). The exact solution is a sphere of radius 0.5√3 centered at (0.5, 0.5, 0.5). The radius of the computed sphere was only 1.3% greater than optimal. I repeated my experiment 10 times and the worst fit was 2.2% greater than optimal. It seems Ritter’s cube experiment result was either unlucky or flawed, and in either case pessimistic.

Sphere example

Next I wanted to generate around 10,000 points uniformly from a sphere of radius 0.5 centered at (0.5, 0.5, 0.5). I did this by generating points from the unit cube and throwing away those that weren’t in the sphere I was after [2]. The sphere takes up about half of the unit cube, so I generated 20,000 points in order to have about 10,000 that landed inside the sphere.

    def inside(c):
        return (c[0]-0.5)**2 + (c[1]-0.5)**2 + (c[2]-0.5)**2 < 1/4

    cube = uniform.rvs(size=(20000, 3)).tolist()
    x = [c for c in cube if inside(c)]
    x[-1] = [0.5, 0.5, 0.0]
    x[-2] = [0.0, 0.5, 0.5]
    x[-3] = [0.5, 0.0, 0.5]
    x[-4] = [0.5, 1.0, 0.5]
    print( bounding_sphere(x) )

This returned a sphere of radius 0.5282 centered at (0.4915, 0.4717, 0.4955). This radius is about 5.6% larger than the optimal radius of 0.5.

I repeated the experiment 10 times, and on average the spheres were 5% larger than optimal. The worst case was 7.6% larger than optimal.

Update

Thanks to Chris Elion for pointing out that the implementation of Ritter’s algorithm that I found does not fully implement the algorithm: it does not update the center as Ritter does in [1].

Chris reran the calculations above and in each case roughly cut the distance between the computed and optimal radius in half.

Related posts

[1] Graphics Gems. Andrew S. Glassner, editor. Academic Press. 1990.

[2] There’s an apocryphal quote attributed to Michelangelo that when asked how he carved David, he said that he started with a block of marble and cut away everything that wasn’t David. In probability this approach is known as acceptance-rejection sampling. It’s simple, but not always the most efficient approach.

6 thoughts on “Approximate minimal bounding sphere

  1. The first started with a sphere “randomly populated with 100,000 points.” Random according to what distribution? I assume uniform. The calculated sphere had a radius 4% larger than optimal. The second example was a cube “randomly populated with 10,000 points.” Again I assume the distribution was uniform.

    This doesn’t solve anything. “Random” means “uniform” by convention, sure. But changing the word leaves you with the same problem you already had — what does the word mean? What do you mean by a random point within a sphere? The closely-analogous problem of “uniformly selecting” a point within a circle is known to be ill-defined. You need to decide what uniform selection means before you can do it.

    https://en.wikipedia.org/wiki/Bertrand_paradox_(probability)

  2. Responding to Michael Watts: I might be being naive, I don’t see that selecting a point on a disc or ball at uniform is ill-defined. It’s straightforward to define a probability distribution such that the probability of picking a point within an area (or volume) is proportional to the size of that area (or volume). You do have to be careful; for example, picking an angle at uniform between 0 and 2 pi and then a distance at uniform between 0 and the radius doesn’t work. But you can show that it doesn’t work, and demonstrate a method that does.

    The Bertrand paradox runs into problems because there’s not a single natural “space of chords” of “uniform density” to sample from.

  3. @Dan: Right, the uniform distribution on a disc or ball, or any bounded and reasonably well-behaved region, is perfectly well-defined. The simplest (although not most efficient) way to generate it is to draw a box around the region and generate uniformly distributed points until you get one inside the region.

    The Bertrand paradox is not really paradoxical at all; it just amounts to the observation that different methods of generating something at random lead to different distributions, and that if one random variable (say, a point inside a circle) is uniformly distributed, another one that is a function of it (say, the distance from that point to the center of the circle) likely is not.

  4. The Bertrand paradox runs into problems because there’s not a single natural “space of chords” of “uniform density” to sample from.

    But there is a clean one-to-one relationship between non-diameter chords and non-origin points within the circle. If you can choose from one of those spaces “at random”, then you can also choose from the other one “at random”, because there’s no difference.

    But we know we cannot choose from the space of non-diameter chords “at random”. The only possible conclusion is that we also cannot choose from the space of non-origin points within a circle “at random”.

    You can define a distribution for which the probability of a point within any region being selected is proportional to the area of the region. But you can’t say that that distribution is “more uniform” than some other distribution.

    The aesthetically cleanest expression of the chord space presented in the wikipedia article is obviously the method of choosing two points on the circumference and drawing the chord between them – that method generates every chord, including each diameter, with equal probability.

    If you’re comfortable with the area-proportional selection of a point within a disc being uniquely “the uniform distribution” over the points within the disc, why aren’t you comfortable with the two-point chord selection model being uniquely “the uniform distribution” over the chords crossing the disc?

  5. @Michael: There’s also a clean one-to-one relationship between points on a unit disc and points in a rectangle with sides 1 and 2pi (use the r and theta of your disc points), but a uniform distribution in one does not look uniform in the other, and I’m fine with that. I am perfectly comfortable with different definitions of uniform distribution for the two cases that depend on their geometry. Probably we are talking past each other at this point but that’s okay. I do agree with you that it is important to be precise about what a uniform distribution actually means when you use the word, but I also think that the “naive” method of defining a uniform distribution on some area/volume is the clear standard. (For example, exercises in probability textbooks will routinely involve uniform distributions over areas/volumes without defining them precisely to your standard every time.)

Leave a Reply

Your email address will not be published.