c++ – Welzl Algorithm to search out the Smallest Bounding Sphere

Your base instances for two, 3, and 4 boundary factors are incorrect.

• For 2 boundary factors, you need the middle to be half the sum of the 2
  boundary factors (ie. their common), not half their distinction. The radius will likely be half the size of their distinction, not 1 / 4. We will regulate your code as follows:

  const glm::vec3 halfSpan = (sos[1] - sos[0]) * 0.5f;
  return Sphere(sos[0] + halfSpan, glm::size(halfSpan));

• For 3 boundary factors, you need the circumcenter of the triangle, not the centroid as you are computing now. Utilizing the formulation given on Wikipedia right here, we are able to write:

$$start{align}
a &= sos[1] – sos[0]
b &= sos[2] – sos[0]
c &= sos[0]

middle &= frac ^2 b – a occasions b + c
finish{align}$$

It’s not regular for the algorithm to attempt to add greater than 4 boundary factors to the set of assist. Since we solely ever add one boundary level at a time, if a selected invocation of the strategy tries so as to add a fifth boundary level, then it will need to have itself been known as with 4, which suggests it could have fallen into the 4-boundary-point base case above and returned a circumsphere of these 4 factors with none try at an additional recursive name.

You do not want to test whether or not the 4 boundary factors are co-planar. If we attempt to add some extent to our set of assist, then which means we already checked and located the minimal sphere containing the remainder of the set doesn’t comprise that time. This may be true of at most three factors on any given aircraft. Of any set of 4 factors on a aircraft, one have to be inside or on the circumcircle outlined by the opposite three, and so already contained in any sphere containing the opposite three.

Additionally, notice that the algorithm is designed to work by eradicating a random level from the set in every recursion. In the event you make this choice non-randomly, in a set order as you are doing right here, then the efficiency traits confirmed in Welzl’s paper could now not apply, and on pathological information you could possibly get runtimes for much longer than anticipated.

Edit: Here is a working C# implementation I wrote as much as perceive the algorithm. I do know you are utilizing C++, however this could assist disambiguate among the particulars.

public static BoundingBall GetBoundingBall(Checklist factors) {
    if (_boundaryPoints == null)
        _boundaryPoints = new Checklist(4);
    else
        _boundaryPoints.Clear();

    return BallWithBounds(factors);
}

static BoundingBall BallWithBounds(Checklist contained) {
    if (contained.Rely == 0 || _boundaryPoints.Rely == 4) {
        swap (_boundaryPoints.Rely) {
            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 final = contained.Rely - 1;
    int removeAt = Random.Vary(0, contained.Rely);

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

    var ball = BallWithBounds(contained);

    if(!ball.Comprises(eliminated)) {   
        _boundaryPoints.Add(eliminated);
        ball = BallWithBounds(contained);        
        _boundaryPoints.RemoveAt(_boundaryPoints.Rely - 1);
    }

    contained.Add(eliminated);
    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 middle = c + Vector3.Cross(A.sqrMagnitude * B - B.sqrMagnitude * A, cross)
                   / (2f * cross.sqrMagnitude);

    float radius = Vector3.Distance(a, middle);

    return new BoundingBall(middle, radius);
}

static BoundingBall TetrahedronCircumSphere(
          Vector3 p1,
          Vector3 p2,
          Vector3 p3,
          Vector3 p4
) { 
    // Assemble a matrix with the vectors as columns 
    // (Xs on one row, Ys on the subsequent... and the final 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 are able to modify it 
    // and nonetheless learn rows from the unique.
    var D = matrix;
    Vector3 middle;

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

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

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

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

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