But calculating a minimum bounding sphere for a point cloud is not a fast operation - it's worse than linear time (which matters for huge meshes that are loaded while the simulator is running) and algorithms usually require robust mathematical operators.
Instead I use a bit of a heuristic - an incremental grow algorithm:
Initialize the sphere to size 0 and center of the first point.
For each point until there are none left.
If the point is inside the bounding sphere, ignore it.
Otherwise, grow the bounding sphere to just barely include this point.This algorithm is pretty good but not optimal. The reason that it's not optimal is that when we grow to encompass a new point P, the opposite point that is "growing" (and not moving) the sphere is the farthest point on the sphere from P. This farthest point may not actually be an input point at all - it could just be an artifact of the fact that we are using the sphere instead of the original data to grow.
(A slower but more comprehensive algorithm would have us go back over the original point data to find the far-side point limits, giving us O(N^2) - plus we'd start to have floating point robustness issues.)
The quality of the bounding sphere has a lot to do with the order in which we grow it - the earliest points establish the size of a sphere, effectively inducing "phantom" points around the temporary (smaller) sphere.
So my first idea was to find the longest axis between any two points in my cloud and insert them first, in the hope that by establishing the long axis we could avoid growing the sphere in false directions. This option probably means fewer grow operations (since more points will be inside the long axis) but the cost of finding the longest axis is O(N^2) so it's not really a speed win.
What surprised me is that using the longest axis vs. the native submit order coming out of the host program, the longest axis was better sometimes but worse in others. The "natural" order of the points seemed to produce pretty good results a lot of the time too.
What I then tried was a randomized approach:
For each of N trials
Shuffle the data
Calculate the bounding sphere
If the sphere is smaller than our previous best
(or this is the first trial)
Save this as our bestNow the quality of these bounding spheres are random (and some are quite poor) but as we increase N, we are more and more likely to randomly find an order that is superior to any pre-determined heuristic. The shuffle takes linear time, and the number of trials is constant (and quite possibly a lot smaller than the number of points).
For N = 256, this produces better results than either the natural or long-axis-based sphere perhaps 90% of the time. And when it doesn't produce the best sphere, it is very close to optimal, usually within a few percent.
This comment has been removed by the author.
ReplyDeleteMaybe I'm misunderstanding the problem (or showing off my poor algorithms skillz)...but can't you do a linear time scan through the points for Max X, Y, and Z and Min X, Y, and Z, and using those values to calculate a bounding sphere?
ReplyDeletefor each point {
if x < min_x: min_x = x
if y < min_y: min_y = y
if z < min_z: min_z = z
if x > max_x: max_x = x
if y > max_y: max_y = y
if z > max_z: max_z = z
}
dx = max_x - min_x
dy = max_y - min_y
dz = max_z = min_z
center = (min_x + dx/2,
min_y + dy/2,
min_z + dz/2
)
bounds = sphere(center, max(dx,dy,dz))
@Jason: That's exactly what I was thinking...
ReplyDelete@Jason: Consider the three points (1, 1), (1.1, 0), and (0, 1.1). Your algorithm would create a bounding sphere with a radius of 1.1, while it should have been sqrt(2). (simplified, you might need more points for this to be true)
ReplyDeleteOn the other hand, Nimrod Megiddo published an algorithm in 1984, which creates a bounding sphere in linear time.
@BioTronic: Aha! Thanks for showing me the error of my ways.
ReplyDelete@BioTronic: Which paper is it?
ReplyDelete@jsnx: That would be "Linear-Time Algorithms for Linear Programming in R³ and Related Problems", available at http://theory.stanford.edu/~megiddo/pdf/lp3.pdf
ReplyDeleteThe centre of the sphere is the average x/y/z of all the points, the radius is from the point furthest from this center - thats O(N) and with a bit more maths can be done in a single pass.
ReplyDelete@baby-rabbit: Your assumption that the centre of the sphere is the average of the points is wrong. If you tried to prove it you'd quickly find counterexamples...
ReplyDeleteMegiddo is linear time but - look at the geometric primitives in each step...it's going to be a slow linear time.
ReplyDeleteWelzl is more useful in practice - see here.
http://www.gamedev.net/reference/programming/features/welzlminsphere/
What I am doing can be thought of as sort of "imperfect in optimization" instead of "imperfect in time".
Basically my algorithm skips the expensive restart steps and uses the incorrect (but close) incremental addition prt of Welzl's algorithm, but by running on multiple permutations, I can get luckier about what points define the sphere.
Check out Timothy Chan's "A simple streaming algorithm for minimum enclosing balls." It is similar to your algorithm, would be just as quick, but it guarantees a 3/2 (radius) approximation to the minimum enclosing ball. The algorithm is dead simple and works fantastically.
ReplyDeletetixxit - the algo in that paper _is_ X-Plane's "cheap" (one-pass) sphere algo. I'm not sure if the algo is original to Chan - we developed it internally but never bothered to do rigorous analysis on its quality.
ReplyDelete(It was the obvious improvement over first-pt-is-center or centroid-is-center.)
The random permuted version simply runs that algo multiple times in random order.