# Welzl Algorithm to find the Smallest Bounding Sphere

I'm implementing Welzl's Algorithm to find the Smallest Bounding Sphere.

In brief, the algorithm works by maintaining:

• a set of points to contain in the sphere (inside its interior, or on its boundary)
• a set of points known to be on the boundary of the sphere (the "set of support")

It works recursively, starting with the full set of points in the contained set and no boundary points in the set of support. At each level, it tries removing a point from the contained set, and recursively finding the minimal bounding sphere of this reduced set (with the current set of support).

• If the removed point is in/on the resulting sphere, then it returns that result up the recursive chain.

• If not, then the removed point must be on the boundary of the minimal bounding sphere, so it adds it to the set of support and recurses again.

Once it boils down to no remaining contained points or 4 boundary points, the recursion bottoms-out and returns an explicitly calculated sphere based on the boundary points found.

I'm not sure of a few details of the implementation, and haven't found on the net the answer to my questions :

• How can I construct a sphere is uniquely defined by four (non co-planar) points?
• Do I have to check in the algorithm if these four points are non co-planar?

• I used the barycenter (centroid) of the boundary points in by set of support when the number such points is three or more. Is that correct ?

• My algorithm may fall into cases with more than four boundary points in the set of support, so I had to add a default case to handle this - is this normal?

Here's my current code:

Sphere ComputeMinimalBoundingShpere(glm::vec3*              points,
std::size_t             numberPoints,
std::vector<glm::vec3>& sos,
std::size_t             numberSos)
{
// Stop the recursion when we have no more points
// Or the number of support is more than three
if (numberPoints == 0 || numberSos >= 4) {
switch (numberSos) {
case 0: return Sphere();
case 1: return Sphere(sos[0], 0.000f);
case 2: {
const glm::vec3 center = (sos[1] - sos[0]) * 0.5f;
return Sphere(center, glm::length(center) * 0.5f);
};
case 3: {
// If we have three points: Look for the barycenter
const glm::vec3& A = sos[0];
const glm::vec3& B = sos[1];
const glm::vec3& C = sos[2];
const float      v = 1.0f / 3.0f;
const float      w = 1.0f / 3.0f;

const glm::vec3 P      = ((1 - v - w) * A) + v * B + w * C;
const float     radius = glm::length(P - A);
};
case 4: {
// If we have four points: Look for the barycenter
const glm::vec3& A = sos[0];
const glm::vec3& B = sos[1];
const glm::vec3& C = sos[2];
const glm::vec3& D = sos[3];
const float      v = 0.25f;
const float      w = 0.25f;
const float      x = 0.25f;

const glm::vec3 P = ((1 - v - w - x) * A) + v * B + w * C + x * D;
const float     radius = glm::length(P - A);
};
default:
assert(false && "We cannot have more than four Set of Support");
}
}

const int    indexPoint = numberPoints - 1;
const Sphere smallestSphere =
ComputeMinimalBoundingShpere(points, numberPoints - 1, sos, numberSos);
if (smallestSphere.contains(points[indexPoint])) { return smallestSphere; }

sos[numberSos] = points[indexPoint];
return ComputeMinimalBoundingShpere(points, numberPoints - 1, sos,
numberSos + 1);
}

• Please clarify your specific problem or add additional details to highlight exactly what you need. As it's currently written, it’s hard to tell exactly what you're asking. – Kromster Aug 12 '18 at 13:53

Your base cases for 2, 3, and 4 boundary points are incorrect.

• For two boundary points, you want the center to be half the sum of the two boundary points (ie. their average), not half their difference. The radius will be half the length of their difference, not a quarter. We can adjust your code as follows:

const glm::vec3 halfSpan = (sos[1] - sos[0]) * 0.5f;
return Sphere(sos[0] + halfSpan, glm::length(halfSpan));
• For three boundary points, you want the circumcenter of the triangle, not the centroid as you're computing now. Using the formulas given on Wikipedia here, we can write:

\begin{align} a &= sos[1] - sos[0]\\ b &= sos[2] - sos[0]\\ c &= sos[0]\\ \\ center &= \frac {\left(\|a\|^2 b - \|b\|^2 a \right) \times \left( a \times b \right)}{2 \|a \times b\|^2} + c \end{align}

It is not normal for the algorithm to try to add more than four boundary points to the set of support. Since we only ever add one boundary point at a time, if a particular invocation of the method tries to add a fifth boundary point, then it must have itself been called with four, which means it would have fallen into the 4-boundary-point base case above and returned a circumsphere of those four points without any attempt at a further recursive call.

You do not need to check whether the 4 boundary points are co-planar. If we try to add a point to our set of support, then that means we already checked and found the minimal sphere containing the rest of the set does not contain that point. This can be true of at most three points on any given plane. Of any set of four points on a plane, one must be inside or on the circumcircle defined by the other three, and so already contained in any sphere containing the other three.

Also, note that the algorithm is designed to work by removing a random point from the set in each recursion. If you make this selection non-randomly, in a fixed order as you're doing here, then the performance characteristics proven in Welzl's paper may no longer apply, and on pathological data you could get runtimes much longer than expected.

Edit: Here's a working C# implementation I wrote up to understand the algorithm. I know you're using C++, but this should help disambiguate some of the details.

public static BoundingBall GetBoundingBall(List<Vector3> points) {
if (_boundaryPoints == null)
_boundaryPoints = new List<Vector3>(4);
else
_boundaryPoints.Clear();

return BallWithBounds(points);
}

static BoundingBall BallWithBounds(List<Vector3> contained) {
if (contained.Count == 0 || _boundaryPoints.Count == 4) {
switch (_boundaryPoints.Count) {
case 0:
return BoundingBall.EMPTY;
case 1:
return new BoundingBall(_boundaryPoints[0], 0);
case 2:
var halfSpan = 0.5f * (_boundaryPoints[1] - _boundaryPoints[0]);
return new BoundingBall(
_boundaryPoints[0] + halfSpan,
halfSpan.magnitude
);
case 3:
return TriangleCircumSphere(
_boundaryPoints[0],
_boundaryPoints[1],
_boundaryPoints[2]
);
case 4:
return TetrahedronCircumSphere(
_boundaryPoints[0],
_boundaryPoints[1],
_boundaryPoints[2],
_boundaryPoints[3]
);
}
}

int last = contained.Count - 1;
int removeAt = Random.Range(0, contained.Count);

Vector3 removed = contained[removeAt];
contained[removeAt] = contained[last];
contained.RemoveAt(last);

var ball = BallWithBounds(contained);

if(!ball.Contains(removed)) {
ball = BallWithBounds(contained);
_boundaryPoints.RemoveAt(_boundaryPoints.Count - 1);
}

return ball;
}

static BoundingBall TriangleCircumSphere(Vector3 a, Vector3 b, Vector3 c) {
Vector3 A = a - c;
Vector3 B = b - c;

Vector3 cross = Vector3.Cross(A, B);

Vector3 center = c + Vector3.Cross(A.sqrMagnitude * B - B.sqrMagnitude * A, cross)
/ (2f * cross.sqrMagnitude);

}

static BoundingBall TetrahedronCircumSphere(
Vector3 p1,
Vector3 p2,
Vector3 p3,
Vector3 p4
) {
// Construct a matrix with the vectors as columns
// (Xs on one row, Ys on the next... and the last row all 1s)
Matrix4x4 matrix = new Matrix4x4(p1, p2, p3, p4);
matrix.SetRow(3, Vector4.one);

float a = matrix.determinant;

// Copy the matrix so we can modify it
// and still read rows from the original.
var D = matrix;
Vector3 center;

Vector4 squares = new Vector4(
p1.sqrMagnitude,
p2.sqrMagnitude,
p3.sqrMagnitude,
p4.sqrMagnitude
);

D.SetRow(0, squares);
center.x =  D.determinant;

D.SetRow(1, matrix.GetRow(0));
center.y = -D.determinant;

D.SetRow(2, matrix.GetRow(1));
center.z =  D.determinant;

center /= 2f * a;
return new BoundingBall(center, Vector3.Distance(p1, center));
}

• Thanks for your answer, it's really valuable. I managed to have the first three spheres but still have an issue with the Tetrahedron. Looking at your code, I don't understand why you have the Get/Set Row shouldn't be Column? As in the Wolfram example? – FloFu Aug 15 '18 at 9:53
• Transposing a matrix doesn't change its determinant, and it was easier to construct the matrix with columns first and then overwrite rows rather than to do the opposite. Same result, just a different intermediate layout for fewer lines of code. – DMGregory Aug 15 '18 at 11:12
• The issue I have on my case is the determinant a is most often equal to 0 (or close to epsilon). Therefore my sphere have an infinite center and radius. – FloFu Aug 15 '18 at 12:03
• I had similar errors when I was constructing the matrix incorrectly. Try calling this method with four very distinct vectors, printing out the entries of the matrix you're using at each step to make sure the layout matches the Wolfram examples (or their transpose) – DMGregory Aug 15 '18 at 12:07
• Not succeed yet :(. Determinant equals zero are due to sos that are equals (maybe an issue with my algorithm). I temporarily break down to just the Sphere From Tetrahedron but it give me invalid results. From the debugger point of view every thing appears to be correct (or I miss something) : pastebin.com/uXPQSv4e – FloFu Aug 15 '18 at 16:31