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
defaultcase 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);
return Sphere(P, radius);
};
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);
return Sphere(P, radius);
};
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);
}
