Hi everyone,

I'm working on implementing Welzl's algorithm for finding the minimal bounding sphere of a set of vertices. I'm getting some pretty huge bounding spheres though, which can't be right. I'm thinking that something must be wrong with how I'm calculating my sphere from either 3 points or from 4 points (possibly both are wrong). As far as I can tell, there's 4 different ways to build a sphere:

sphere from 1 point: Center is point, radius is 0

sphere from 2 points: Center is midpoint, radius is half distance

sphere from 3 points: I use this method:


Sphere::Sphere( Vector& p1, Vector& p2, Vector& p3 )
{
Vector v23 = p2 - p3;
Vector v13 = p1 - p3;
Vector v12 = p1 - p2;
Vector v21 = p2 - p1;
Vector v31 = p3 - p1;
Vector v32 = p3 - p2;

float denominator = ( 2 * ( v12.CrossProduct(v23).Length() * v12.CrossProduct(v23).Length() ) );
float alpha = ( (v23.Length() * v23.Length()) * (v12.DotProduct(v13)) ) / denominator;
float beta = ( (v13.Length() * v13.Length()) * (v21.DotProduct(v23)) ) / denominator;
float gamma = ( (v12.Length() * v12.Length()) * (v31.DotProduct(v32)) ) / denominator;

center = alpha * p1 + beta * p2 + gamma * p3;

radius = ( v12.Length() * v23.Length() * v31.Length() ) / ( 2 * v12.CrossProduct(v23).Length() );
}


I've tried several methods for the 4 point sphere, here's the one that I've settled on:
4 point sphere method

Now when I expand a sphere, I take the support set (all points currently used to define the sphere, always less than 5) and the new point and test all possible permutations of the support set and the new point to find the minimal sphere. For instance, if my support set has 4 vertices then I first try all two point combinations of the support set and the new point, then all 3 point combinations and finally all 4 point combinations. I then replace/remove any necessary points from the support set and add the newpoint to the support set. I'm 99% sure that the logical flow of my algorithm is right, I believe I just have some grievous math error somewhere in my Sphere calculation routine.

Are there any better ways to calculate a minimal bounding sphere from 3 or 4 points?

Thank you.

Edit: Also, on the off chance that anyone needs to see my sphere calculation functions in order to help:

Full Sphere Calculation Functions

Let your points be x1..4. Then the center of your sphere is the c which satisfies the following 3 equations:
|x2 - c|^2 = |x1 - c|^2
|x3 - c|^2 = |x1 - c|^2
|x3 - c|^2 = |x1 - c|^2

Multiplying out

|x2|^2 - 2(x2 `dot` c) + |c|^2 = |x1|^2 - 2(x1 `dot` c) + |c|^2

and similarly for the other 2 equations

Rearranging this gives

2(x1 - x2) `dot` c = |x1|^2 - |x2|^2
2(x1 - x3) `dot` c = |x1|^2 - |x3|^2
2(x1 - x4) `dot` c = |x1|^2 - |x4|^2

Which is a simple system of 3 linear equations in 3 unknowns, which can be solved using Gaussian elimination for example.

Once you have the center, the radius is trivial to calculate.

quasar3d: Do you realize that the sphere that touches four points might be much larger than the smallest sphere that contains those four points? For instance, imagine a sphere of radius 6,378 Km and pick four points on the surface that are near each other (say, within 1m of each other). The sphere that passes through all four of them will be the original Earth-size sphere, but you can certainly find much smaller spheres that surround all four points.

As for the OP, if you understand the method you are using and you have a specific example that fails, you should follow the computations with a debugger and do them on your own in parallel, to see where the code is going wrong.

@alvero, I know that:), but if this sphere isn't minimal, then the minimal sphere is the sphere trough a subset of the points. I assume that's why lukabratzi mentions those other cases as well.

lukabratzi: Of course it's easy enough to test if your 4 vertices actually lie on your sphere. If this isn't true, than it's pretty clear the math parts are wrong, and otherwise it must be your algorithm.

Rearranging this gives

2(x1 - x2) `dot` c = |x1|^2 + |x2|^2
2(x1 - x3) `dot` c = |x1|^2 + |x3|^2
2(x1 - x4) `dot` c = |x1|^2 + |x4|^2

Which is a simple system of 3 linear equations in 3 unknowns, which can be solved using Gaussian elimination for example.

Once you have the center, the radius is trivial to calculate.

Thanks for the help. I do have a question though: how can I take the magnitude squared of x1, x2, etc? Technically, they're points not vectors right?

[quote name='quasar3d' timestamp='1320417294' post='4880474']
Rearranging this gives

2(x1 - x2) `dot` c = |x1|^2 + |x2|^2
2(x1 - x3) `dot` c = |x1|^2 + |x3|^2
2(x1 - x4) `dot` c = |x1|^2 + |x4|^2

Which is a simple system of 3 linear equations in 3 unknowns, which can be solved using Gaussian elimination for example.

Once you have the center, the radius is trivial to calculate.

Thanks for the help. I do have a question though: how can I take the magnitude squared of x1, x2, etc? Technically, they're points not vectors right?
[/quote]

Just treat them as vectors. In the original equation, you are taking differences between points, which gives a vector, and at that point, treating your points as vectors clearly makes no difference. After that it's just algebra, and pretty much loses its geometric meaning anyway.

I just spotted a mistake:

2(x1 - x2) `dot` c = |x1|^2 + |x2|^2
2(x1 - x3) `dot` c = |x1|^2 + |x3|^2
2(x1 - x4) `dot` c = |x1|^2 + |x4|^2

should be


2(x1 - x2) `dot` c = |x1|^2 - |x2|^2
2(x1 - x3) `dot` c = |x1|^2 - |x3|^2
2(x1 - x4) `dot` c = |x1|^2 - |x4|^2

I just spotted a mistake:

2(x1 - x2) `dot` c = |x1|^2 + |x2|^2
2(x1 - x3) `dot` c = |x1|^2 + |x3|^2
2(x1 - x4) `dot` c = |x1|^2 + |x4|^2

should be


2(x1 - x2) `dot` c = |x1|^2 - |x2|^2
2(x1 - x3) `dot` c = |x1|^2 - |x3|^2
2(x1 - x4) `dot` c = |x1|^2 - |x4|^2

Thanks for all the help. I hate asking you this as it's something that I know I should know, but can you explain a bit more about how to put the above equations into matrix form? The dot product of the vectors with the center is throwing me.

Again, thank you for the help

[quote name='quasar3d' timestamp='1320425130' post='4880514']
I just spotted a mistake:

2(x1 - x2) `dot` c = |x1|^2 + |x2|^2
2(x1 - x3) `dot` c = |x1|^2 + |x3|^2
2(x1 - x4) `dot` c = |x1|^2 + |x4|^2

should be


2(x1 - x2) `dot` c = |x1|^2 - |x2|^2
2(x1 - x3) `dot` c = |x1|^2 - |x3|^2
2(x1 - x4) `dot` c = |x1|^2 - |x4|^2

Thanks for all the help. I hate asking you this as it's something that I know I should know, but can you explain a bit more about how to put the above equations into matrix form? The dot product of the vectors with the center is throwing me.

Again, thank you for the help
[/quote]

A multiplication of a 3x3 matrix with a column 3-vector is simply short-hand for the computation of the dot product of each of the three rows of the matrix with the column. So just build a 3x3 matrix whose rows are the vectors 2(x1-x2), 2(x1-x3) and 2(x1-x4). Is that what you needed?