Bouncing Bubble Algorithm for smallest enclosing sphere

Question 1

I realize this post is a year old, but I implemented the bouncing bubble algorithm in C++ and thought maybe some people would find it useful. This is based on the paper by Tian Bo, which I purchased for $18.

BoundSphere calculateBoundSphere(vector<Vertex> vertices){
    Vector3D center = vertices[0].position;
    float radius = 0.0001f;
    Vector3D pos, diff;
    float len, alpha, alphaSq;
    vector<Vertex>::iterator it;

    for (int i = 0; i < 2; i++){
        for (it = vertices.begin(); it != vertices.end(); it++){
            pos = it->position;
            diff = pos - center;
            len = diff.length();
            if (len > radius){
                alpha = len / radius;
                alphaSq = alpha * alpha;
                radius = 0.5f * (alpha + 1 / alpha) * radius;
                center = 0.5f * ((1 + 1 / alphaSq) * center + (1 - 1 / alphaSq) * pos);
            }
        }
    }

    for (it = vertices.begin(); it != vertices.end(); it++){
        pos = it->position;
        diff = pos - center;
        len = diff.length();
        if (len > radius){
            radius = (radius + len) / 2.0f;
            center = center + ((len - radius) / len * diff);
        }
    }

    return BoundSphere(center, radius);
}

Question 2

Here is the implementation you asked for. For analysis part, you can get the idea reading the following.

These anaylsis are for circle(2d) but the extension to sphere(3d) should be easy. I am sure the analysis will help.

An O(n2)-time Algorithm

Draw a circle at center, c, such that points of given set lie within the circle. Clearly, this circle can be made smaller (or else we have the solution).
Make a circle smaller by finding the point A farthest from the center of circler, and drawing a new circle with the same center and passing through the point A. These two steps produce a smaller enclosing circle. The reason that the new circle is smaller is that this new circle still contains all the points of given set, except now it passes through farthest point, x, rather than enclosing it.
If the circle passes through 2 or more points, proceed to step 4. Otherwise, make the circle smaller by moving the center towards point A, until the circle makes contact with another point B from the set.
At this stage, we a circle, C, that passes through two or more points of a given set. If the circle contains an interval (point-free interval) of arc greater than half the circle's circumference on which no points lie, the circle can be made smaller. Let D and E be the points at the ends of this point-free interval. While keeping D and E on the circle's boundary, reduce the diameter of the circle until we have either case (a) or case (b).
- Case (a) The diameter is the distance DE.
  We are done!
- Case (b) The circle C touches another point from the set, F.
  Check whether there exits point-free arc intervals of length more than half the circumference of C.
  IF
  no such point-free arc intervals exit THEN We are done!
  Else
  Goto step 4. In this case, three points must lie on an arc less than half the circumference in length. We repeat step 4 on the outer two of the three points on the arc.

Analysis

This algorithm is straight forward, but expensive to implement. Step 1, 2 and 3 require linear time in the number of points in the given set. In step 4 above, it takes linear time to find each new point F. However, finding the point F does not guarantee the termination of the algorithm. Step 4 must be repeated until the circle contains no point-free interval longer than half its circumference. In the worst case, this requires (n-2) iterations of step 4, implying that the total time spent in step 4 could be on the order of the square of the size of the point set.

Hence, this algorithm is O(n2).

O(n) Algorithm

Nimrod Megiddo proposes the algorithm and it uses the prune-and-search techniques for linear programming to find the minimal enclosing circle in O(n) time.

Prune-and-Search

The essence of Megiddo's algorithm is prune-and-search. In a prune-and-search algorithm, a linear amount of work is done at each step to reduce the input size by a constant fraction f. If this can be achieved, then the total amount of work done will reduce to O(n)*(1 + (1-f) + (1-f)2 + ...). In this formula, the infinite series is geometric and sums to a constant value, and so the total running time is O(n).

For example, suppose that by inspecting our set of n input elements we can discard 1/4 of them as irrelevant to the solution. By repeatedly applying this inspection to the remaining elements, we can reduce the input to a size which is trivial to solve, say n≤3. The total time taken to achieve this reduction will be proportional to (n + 3n/4 + 9n/16 + ...). It is easy to show that this series approaches, and never exceeds, a limit of 4n. Therefore, the total running time is proportional to n, as required.

The idea of using the geometric series to reduce an algorithm to linear time predates Megiddo's work; in particular, it had previously been used to develop O(n) median-finding algorithms. However, he was the first to apply it to a number of fundamental problems in computational geometry.

Megiddo's Linear-Time Algorithm

To find the minimal enclosing circle (MEC) of a set of points, Megiddo's algorithm discards at least n/16 points at each (linear-time) iteration. That is, given a set S of n points, the algorithm identifies n/16 points that can be removed from S without affecting the MEC of S. This procedure can be repeatedly applied until some trivial base case is reached (such as n=3), with the total running time proportional to (n + 15n/16 + 225n/256 + ...) = 16n.

In order to find n/16 points to discard, a great deal of cleverness is required. The algorithm makes heavy use of two subroutines:

median(S, >): Takes a set S of elements and a metric for comparing them pairwise, and returns the median of the elements.
MEC-center(S, L): Takes a set S of points and a line L, and determines which side of L the center of the MEC of S lies on.

As mentioned above, median() predates Megiddo's work, whereas the algorithm described here as MEC-center() was presented as part of his 1983 paper. To explore these procedures in detail would go beyond the scope of this outline, but each uses prune-and-search to run in linear time. The algorithm used by MEC-center() amounts to a simplified version of the algorithm as a whole.

Given these primitives, the algorithm for discarding n/16 input points runs as follows:

Arbitrarily pair up the n points in S to get n/2 pairs.
Construct a bisecting line for each pair of points, to get n/2 bisectors in total.
Call median() to find the bisector with median slope. Call this slope mmid.
Pair up each bisector of slope ≥ mmid with another of slope < mmid, to get n/4 intersection points. Call the set of these points I.
Call median() to find the point in I with median y-value. Call this y-value ymid.
Call MEC-center() to find which side of the line y=ymid the MEC-center C lies on. (Without loss of generality, suppose it lies above.)
Let I' be the subset of points of I whose y-values are less than ymid. (I' contains n/8 points.)
Find a line L with slope mmid such that half the points in I' lie to L's left, half to its right.
Call MEC-center() on L. Without loss of generality, suppose C lies on L's right.
Let I'' be the subset of I' whose points lie to the left of L. (I'' contains n/16 points.)

We can now discard one point in S for each of the n/16 intersection points in I''. The reasoning runs as follows. After two calls to MEC-center(), we have found that the MEC-center C must lie above ymid and to the right of L, whereas any point in I'' is below ymid and to the left of L.

Each point in I'' is at the meeting point of two bisector lines. One of these bisectors must have slope ≥ mmid, and therefore must never pass through the quadrant where we know C to lie. Call this bisector B. Now, we know which side of B C lies on, so of the two points whose bisector is B, let PC be the one that lies on the same side as C, and let the other be PX.

It is easy to show that PC must be nearer to C than PX. It follows that PC cannot lie on the minimal enclosing circle, and thus we can safely discard a point PC for each of the n/16 intersection points in I''.

We have not discussed here how this algorithm can be made to deal with degenerate input (parallel bisectors, collinear points, etc), but it turns out that we get the same performance guarantees for such cases. The fact of the matter is for degenerate input the algorithm is able to discard more than n/16 points. In short, Megiddo's algorithm guarantees to prune at least n/16 points in each iteration independent to input.

Therefore, by the argument based on the geometric series, Megiddo's algorithm computes the minimal enclosing circle in linear time.

Refer to this paper for more on the O(n) Algorithm

Question 3

UPDATE: Upon closer review of the animation and the Wikipedia article, I've decided my algorithm is not Bouncing Bubble. From the animation it is clear that Bouncing Bubble moves the bounding sphere, while my algorithm only grows it in one direction at a time. Also, my understanding is that Bouncing Bubble uses an approximation and may not contain all the points; however, there seems to be a tolerance parameter that you can control.

I think my algorithm still has some nice properties:

It's O(n).
It's incremental.
It's guaranteed to contain all the points.
It's dead simple. The only thing simpler would be to find the axis-aligned bounding box, and then put a sphere around that, but that's surely going to be a larger volume than this algorithm produces.

I came up with a O(n) incremental algorithm on my own, which from the description on Wikipedia I believe must be the Bouncing Bubble algorithm. I have no idea how close it comes to the minimally bounding circle. First I'll try to explain my reasoning behind the algorithm. It may help to draw this on paper.

Imagine you have a bounding circle with some center C and radius r. You want to extend the circle to contain a new point P. You compute the distance d from C to P and find that it is greater than r, so it must be outside the circle. Now if you cast a ray from P through C, it exits the back of the circle at A.

Now imagine you pin the circle down at A and then expand it along the ray back toward P. It will continue to contain all the previous points, and now it contains P. The diameter of this new circle is r + d, thus its radius is r' = (r + d) / 2. Walk along the ray from P to A by the distance of the new radius, and you're at the center of the new circle.

Here's a sketch of the algorithm:

Initialize the bounding sphere to be centered on the first point C with a radius r of 0.
For each remaining point P:
1. If the distance d from P to C is <= r, skip it.
2. If d > r:
  1. The new radius r' = (r + d) / 2.
  2. The new center C' = P + r' (C - P) / d.

And here's some Objective-C code I wrote just today:

- (void)updateBoundingVolume {
    self.boundingSphereCenter = GLKVector3Make(0, 0, 0);
    self.boundingSphereRadius = 0;

    GLfloat* vertex = self.vertices;
    GLfloat* endV = vertex + 3 * self.vertNum;

    if (vertex < endV) {
        // initialize to zero radius sphere centered on the first point
        self.boundingSphereCenter = GLKVector3Make(vertex[0], vertex[1], vertex[2]);
        vertex += 3;
    }

    while (vertex < endV) {
        GLKVector3 p = GLKVector3Make(vertex[0], vertex[1], vertex[2]);
        vertex += 3;

        float d = GLKVector3Distance(self.boundingSphereCenter, p);
        if (d <= self.boundingSphereRadius) {
            // already inside the sphere
            continue;
        }

        // length of line from other side of sphere through the center to the new point
        // divide by 2 for radius
        self.boundingSphereRadius = (self.boundingSphereRadius + d) / 2;

        // unit vector from new point to old center
        GLKVector3 u = GLKVector3DivideScalar(GLKVector3Subtract(self.boundingSphereCenter, p), d);
        // step back from the new point by the new radius to get the new center
        self.boundingSphereCenter = GLKVector3Add(p, GLKVector3MultiplyScalar(u, self.boundingSphereRadius));
    }
}