Reputation: 10655
I'm looking for an algorithm which someone has access to that will compute the smallest bounding sphere that encloses a set of other bounding spheres. I have thought about this for a while and have come up with some initial solutions, but I don't believe these are necessarily the most accurate or the least computationally expensive (fastest).
My first solution is the simplest naive one, which is to average out the sphere centers to get the center point, and then compute the maximum distance from the calculated center to each sphere's center plus its radius, as the radius. So pseudo code goes like:
function containing_sphere_1(spheres)
center = sum(spheres.center) / count(spheres)
radius = max(distance(center, spheres.center) + radius)
return Sphere(center, radius)
end
However I get the feeling that it isn't that computationally cheap, nor is it quite accurate since the resulting sphere could be quite larger than it needs to be.
My second thought is to use an iterative algorithm to compute the minimal bounding sphere. It is computed by successively testing another sphere, if the tested sphere is inside the bounds, then nothing is done, otherwise a new bounding sphere is computed from the two spheres available. The new bounding sphere has a center that is half way between the vector between the two centers if it was extended to the spheres surfaces, and the radius is the half the length of that line (from the new center to either sphere's surface).
function containing_sphere_2(spheres)
bounds = first(spheres)
for each sphere in spheres
if bounds does not contain sphere
line = vector(bounds.center, sphere.center)
extend(line, bounds.radius)
extend(line, sphere.radius)
center = midpoint(line)
radius = length(line) / 2
bounds = Sphere(center, radius)
end
end
return bounds
end
Initially I thought that this would be the way to go, since it is iterative and seems fairly logically consistent, however after doing some reading, most notably the article "Smallest enclosing disks (balls and ellipsoids)" by Emo Welzl I'm not not so sure.
As I understand it the basis of this algorithm is that the minimum bounding sphere over a set of points in 3 dimensions can be determined by at most 4 points (which are on the surface of the enclosing sphere). So the algorithm takes an iterative approach by selecting 4 points, and then testing other points to see if they're inside or not, if they aren't a new bounding sphere is constructed featuring the new point.
Now the algorithm deals strictly with points, but I think it can be applied to deal with spheres, the main complication being accounding for the radius when constructing the enclosing sphere.
So what is the 'best', as in least computationally expensive, algorithm that creates a minimal bounding sphere for a set of given spheres?
Is one of these I have described here the answer? Some pseudo code or the algorithm would be great.
Upvotes: 15
Views: 19305
Reputation: 399
This is almost certainly not a computationally-optimal approach*, but it could be the code-optimal approach (i.e. lowest effort to write), depending on what routines you have available, which may be useful to some.
Start with your first approach: find the average of the points (perhaps weighted by the radii), and use that as the starting point for a minimization algorithm like scipy.optimize.minimize, minimizing the sum of the radius plus the distance between the current position (center of your bounding sphere) and the other sphere. In python:
import numpy as np
from scipy.optimize import minimize
# spheres given by centers = np.array([[x0,y0,z0],...])
# and radii = np.array([r0,...])
def dist(pos):
return np.min(radii + np.linalg.norm(centers-pos[np.newaxis,:],axis=1))
result = minimize(dist,x0=np.average(centers,axis=0,weights=radii))
bounding_center = result.x
bounding_radius = result.fun
* Time complexity is O((1+F)*N) for N spheres and F calls to the function by the minimization routine.
Upvotes: 0
Reputation: 411
Here's a fast, near optimal approach based on Ritter's algorithm https://en.wikipedia.org/wiki/Bounding_sphere :
For each sphere, find its min/max x/y/z points. Throw these 6 points into a bucket. When you've done all N spheres, you'll have a bucket of 6N points. Find a bounding sphere for these using any of the known algorithms.
The bounding sphere you get will very likely be a little too small, regardless of algorithm. You could then do the 2nd pass of Ritter's method, but using the backsides of the spheres as the points to test. 'Backside' meaning pt on sphere farthest from current bnd sphere's center. If a sphere's backside pt is outside current bnd sphere, grow bnd sphere to include it.
In addition to the 6 extrema pts, you could initially include the 8 corners of the inscribed cube:
( [+/-]kR, [+/-]kR, [+/-]kR ), where k=sqrt(3)/3. This gives 14 pts that are fairly well-distributed, geodesically.
Upvotes: 5
Reputation: 8823
The step from enclosing points to enclosing spheres is non-trivial, as the discussion of Welzl's algorithm (which works to enclose points) in K. Fischer's thesis explains, "Smallest enclosing ball of balls". See Sec. 5.1.
Chapter 4 presents the "enclosing points" material, then Chapter 5 presents "enclosing spheres".
The algorithm described in Fisher's thesis has been implemented in the CGAL package since release 3.0, if you are just looking for an implementation.
Upvotes: 11