Minimum bounding circle (sphere) for a triangle (tetrahedron)

In Chapter 4 of RTCD I talked about, amongst other things, how to compute the minimum bounding sphere for a set of points. The algorithm I described (Welzl’s algorithm, in Section 4.3.5) requires a pair of functions to return the bounding sphere for three as well as four points. As I didn’t cover these primitive functions in the book, I thought I’d here explain how to compute both.

I’ll start with the bounding circle (or sphere) for a triangle. Let ABC be the triangle. What we want to find is a point P in the plane of ABC that is the center of the smallest circle K bounding ABC, as well as the radius r of this circle. There are two cases to consider:

  1. Two triangle vertices lie on K.
  2. All three triangle vertices lie on K.

The first case has one vertex inside K and two vertices on K. The two vertices on K must lie on the diagonal of K or K would not be minimal, which means the diagonal must correspond to the longest side of ABC.

The second case corresponds to K being the circumcircle of ABC, with P being the circumcenter. (The circumcenter is the unique point equidistant to all three vertices.) Here, again, we have two cases:

  1. P lies inside ABC, which occurs iff ABC is acute.
  2. P lies outside (or on) ABC, which occurs iff ABC is an obtuse (or right) triangle.

P can be expressed as lying in the plane of ABC as:

P = A+s*(B-A)+t*(C-A)

where s and t (and r = 1 - s - t) are the barycentric coordinates of P with respect to ABC.

I won’t attempt to show it here, but it turns out that P lying outside ABC means the circumcircle is not the smallest bounding circle and instead tells us that the smallest bounding sphere is determined by the longest side of ABC. Furthermore, it turns out that the side of ABC that P lies outside of, is the longest side!

What this means is that we can distinguish between the first two cases (two or thee vertices on the boundary of the minimum bounding circle) by simply computing s and t. We have three vertices on the boundary if s > 0, t > 0 and 1 - s - t > 0, and then P is given by P = A+s*(B-A)+t*(C-A). If instead one of the terms is negative (and only one term will be negative) we know which edge that P lies outside of, and therefore also which two triangle vertices lie on the diameter of K!

So, now to computing s and t as they’re the key to determine the minimum bounding circle. The condition that P is equidistant to A, B, and C is expressed as:

Dot(P-B,P-B) = Dot(P-A,P-A)
Dot(P-C,P-C) = Dot(P-A,P-A)

This expresses that the squared distance PB, Dot(P-B, P-B), and the squared distance PC, Dot(P-C, P-C), both must be be the same as the squared distance PA, Dot(P-A, P-A). Inserting the expression for P and collecting terms result in:

2*Dot(A-B,B-A)*s + 2*Dot(A-B,C-A)*t + Dot(A-B,A-B) = 0
2*Dot(A-C,B-A)*s + 2*Dot(A-C,C-A)*t + Dot(A-C,A-C) = 0

Rearranging the expressions to obtain fewer distinct subterms gives:

2*Dot(A-B,B-A)*s + 2*Dot(A-B,C-A)*t = -Dot(A-B,A-B)
2*Dot(A-C,B-A)*s + 2*Dot(A-C,C-A)*t = -Dot(A-C,A-C)
 
-[2*Dot(A-B,B-A)*s + 2*Dot(A-B,C-A)*t] = Dot(A-B,A-B)
-[2*Dot(A-C,B-A)*s + 2*Dot(A-C,C-A)*t] = Dot(A-C,A-C)
 
-2*Dot(A-B,B-A)*s + -2*Dot(A-B,C-A)*t = Dot(A-B,A-B)
-2*Dot(A-C,B-A)*s + -2*Dot(A-C,C-A)*t = Dot(A-C,A-C)
 
2*Dot(B-A,B-A)*s + 2*Dot(B-A,C-A)*t = Dot(B-A,B-A)
2*Dot(C-A,B-A)*s + 2*Dot(C-A,C-A)*t = Dot(C-A,C-A)
 
Dot(B-A,B-A)*s + Dot(B-A,C-A)*t = (1/2)*Dot(B-A,B-A)
Dot(C-A,B-A)*s + Dot(C-A,C-A)*t = (1/2)*Dot(C-A,C-A)

These two equations can be solved using Cramer’s rule to give:

dotABAB = Dot(B-A,B-A)
dotABAC = Dot(B-A,C-A)
dotACAC = Dot(C-A,C-A)
d = dotABAB*dotACAC - dotABAC*dotABAC
s = (1/2)*(dotABAB*dotACAC - dotACAC*dotABAC) / d
t = (1/2)*(dotACAC*dotABAB - dotABAB*dotABAC) / d

Here we note that the scalar d could become zero, in which case A, B, C lie on a line. In this case we obtain the circle by computing the axis-aligned bounding box (AABB) for the three points and pick the box center as the circle center, and the distance from this center to any box corner as the circle radius.

Taking the math for the barycentric coordinates for P and the knowledge that P lying outside the triangle means we need to pick the midpoint on the edge we’re outside, we can now turn everything into code:

void MinimumBoundingCircle(Circle &circle, Point a, Point b, Point c) {
    float dotABAB = Dot(b - a, b - a);
    float dotABAC = Dot(b - a, c - a);
    float dotACAC = Dot(c - a, c - a);
    float d = 2.0f*(dotABAB*dotACAC - dotABAC*dotABAC);
    Point referencePt = a;
    if (Abs(d) <= EPSILON) {
        // a, b, and c lie on a line. Circle center is center of AABB of the
        // points, and radius is distance from circle center to AABB corner
        AABB bbox = ComputeAABB(a, b, c);
        circle.c = 0.5f * (bbox.min + bbox.max);
        referencePt = bbox.min;
    } else {
        float s = (dotABAB*dotACAC - dotACAC*dotABAC) / d;
        float t = (dotACAC*dotABAB - dotABAB*dotABAC) / d;
        // s controls height over AC, t over AB, (1-s-t) over BC
        if (s <= 0.0f) {
            circle.c = 0.5f * (a + c);
        } else if (t <= 0.0f) {
            circle.c = 0.5f * (a + b);
        } else if (s + t >= 1.0f) {
            circle.c = 0.5f * (b + c);
            referencePt = b;
        } else circle.c = a + s*(b - a) + t*(c - a);
    }
    circle.r = Sqrt(Dot(circle.c - referencePt, circle.c - referencePt));
}

This approach directly generalizes to computing the minimum bounding sphere for a tetrahedron ABCD. I won’t go through the full derivation again as it’s straightforward given the above information. The only difference from the above is that we instead start with these expressions

P = A+s*(B-A)+t*(C-A)+u*(D-A)
Dot(P-B,P-B) = Dot(P-A,P-A)
Dot(P-C,P-C) = Dot(P-A,P-A)
Dot(P-D,P-D) = Dot(P-A,P-A)

but otherwise do everything exactly the same as before.

Similar Posts:

Share and Enjoy:
  • Digg
  • del.icio.us
  • Facebook
  • Reddit
  • Slashdot
  • StumbleUpon
  • Technorati
  • LinkedIn

6 Comments »

  1. Fred said,

    August 2, 2007 @ 2:45 am

    Hi Christer,

    I don’t wish to be pedantic but …

    In your code example the line:
    float d = 2.0f*dotABAB*dotACAC - dotABAC*dotABAC;

    has a pair of brackets missing:
    float d = 2.0f * ( dotABAB*dotACAC - dotABAC*dotABAC ) ;

    Best regards,

    Steve.
    p.s. Nice blog - keep the articles rolling.

  2. christer said,

    August 2, 2007 @ 8:17 am

    Oops! Right you are! Thanks for the correction and the feedback Steve, the bug’s been fixed now.

  3. ash said,

    September 5, 2007 @ 9:46 pm

    Thanks for posting this. I’m rewriting a ray tracer I built a while ago, and I really wanted to make the collision detection for the tetrahedrons more efficient.

  4. tsuraan said,

    October 8, 2008 @ 6:29 pm

    In the case where all three vertices of the triangle lie on a line, the bounding circle isn’t infinite, it’s just the circle whose diameter is the largest of |AB|, |AC|, or |BC|, centered on the center of the largest line segment.

    Nice write-up other than that :-)

  5. christer said,

    October 8, 2008 @ 7:40 pm

    Thank you for the correction tsuraan; you are, of course, absolutely right! I believe the easiest (and most robust) way of obtaining the circle for this boundary case is by computing the AABB of the points and return the circle as being the center of this AABB and the radius being the distance between the AABB center and any AABB corner.

    I have corrected the text and the code snippet accordingly.

  6. huang said,

    September 17, 2009 @ 5:10 am

    Thanks for posting this.

RSS feed for comments on this post · TrackBack URI

Leave a Comment

You must be logged in to post a comment.