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& 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);
}
Answer
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 points) {
if (_boundaryPoints == null)
_boundaryPoints = new List(4);
else
_boundaryPoints.Clear();
return BallWithBounds(points);
}
static BoundingBall BallWithBounds(List 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)) {
_boundaryPoints.Add(removed);
ball = BallWithBounds(contained);
_boundaryPoints.RemoveAt(_boundaryPoints.Count - 1);
}
contained.Add(removed);
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);
float radius = Vector3.Distance(a, center);
return new BoundingBall(center, radius);
}
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));
}
No comments:
Post a Comment