Calculate circumsphere of tetrahedron

954 Views Asked by At

Looking to calculate the circumsphere of a tetrahedron for my Delaunay triangulation. Following the explanation I found here (bottom half of page), I have the following code:

         public static Vector3 Circumcentre(Vector3 a, Vector3 b, Vector3 c, Vector3 d)
         {
             // Calculate plane for edge AB
             Vector3 planeABNormal = b - a;                      
             Vector3 planeABPoint = Vector3.Lerp(a, b, 0.5f);

             // Calculate plane for edge AC
             Vector3 planeACNormal = c - a;
             Vector3 planeACPoint = Vector3.Lerp(a, c, 0.5f);

             // Calculate plane for edge CD
             Vector3 planeCDNormal = d - c;
             Vector3 planeCDPoint = Vector3.Lerp(c, d, 0.5f);

             // Calculate line that is the plane-plane intersection between AB and AC
             Vector3 linePoint;
             Vector3 lineDirection;

             // Taken from: http://wiki.unity3d.com/index.php/3d_Math_functions
             PlanePlaneIntersection(out linePoint, out lineDirection, planeABNormal, planeABPoint, planeACNormal, planeACPoint);

             // Calculate the point that is the plane-line intersection between the above line and CD
             Vector3 intersection;

             // Taken from: http://wiki.unity3d.com/index.php/3d_Math_functions
             LinePlaneIntersection(out intersection, linePoint, lineDirection, planeCDNormal, planeCDPoint);

             return intersection;
         }

The problem is, the circumsphere only touches 3 of the 4 points of the tetrahedron. What have I done wrong?

EDIT: Example set (rounded to 1dp):

Input:

a = (0.0, 3.4, 7.3), b = (-1.0, -1.7, -1.8), c = (1.0, -4.5, -1.0), d = (0.0, 1.4, -10.6)

Computed values:

planeABNormal = (-1.0, -5.1, -9.1), planeABPoint = (-0.5, 0.9, 2.8)

planeACNormal = (1.0, -7.9, -8.3), planeACPoint = (0.5, -0.6, 3.2)

planeCDNormal = (-1.0, 5.9, -9.6), planeCDPoint = (0.5, -1.6, -5.8)

linePoint = (1.0, -1.0, 3.6), lineDirection = (-29.6, -17.4, 13.0)

intersection = (1.3, -0.8, 3.5) -> this should be the center of the circumsphere

Circumradius: I calculate the circumradius by getting the distance from the circumcentre to point a.

With this input set / result, when I draw the circumsphere points a, b, and c, lie on the circumference of the sphere but point d sticks way out the side of the sphere

1

There are 1 best solutions below

0
On BEST ANSWER

May have figured it out (tests look good so far):

    public static Vector3 Circumcentre(Vector3 a, Vector3 b, Vector3 c, Vector3 d)
    {
        // Calculate plane for edge AB
        Vector3 planeABNormal = b - a;
        Vector3 planeABPoint = Vector3.Lerp(a, b, 0.5f);

        // Calculate plane for edge AC
        Vector3 planeACNormal = c - a;
        Vector3 planeACPoint = Vector3.Lerp(a, c, 0.5f);

        // Calculate plane for edge BD
        Vector3 planeBDNormal = d - b;
        Vector3 planeBDPoint = Vector3.Lerp(b, d, 0.5f);

        // Calculate plane for edge CD
        Vector3 planeCDNormal = d - c;
        Vector3 planeCDPoint = Vector3.Lerp(c, d, 0.5f);

        // Calculate line that is the plane-plane intersection between AB and AC
        Vector3 linePoint1;
        Vector3 lineDirection1;

        // Taken from: http://wiki.unity3d.com/index.php/3d_Math_functions
        PlanePlaneIntersection(out linePoint1, out lineDirection1, planeABNormal, planeABPoint, planeACNormal, planeACPoint);

        Vector3 linePoint2;
        Vector3 lineDirection2;

        PlanePlaneIntersection(out linePoint2, out lineDirection2, planeBDNormal, planeBDPoint, planeCDNormal, planeCDPoint);

        // Calculate the point that is the plane-line intersection between the above line and CD
        Vector3 intersection;

        // Floating point inaccuracy often causes these two lines to not intersect, in that case get the two closest points on each line 
        // and average them
        // Taken from: http://wiki.unity3d.com/index.php/3d_Math_functions
        if (!LineLineIntersection(out intersection, linePoint1, lineDirection1, linePoint2, lineDirection2))
        {
            Vector3 closestLine1;
            Vector3 closestLine2;

            // Taken from: http://wiki.unity3d.com/index.php/3d_Math_functions
            ClosestPointsOnTwoLines(out closestLine1, out closestLine2, linePoint1, lineDirection1, linePoint2, lineDirection2);

            // Intersection is halfway between the closest two points on lines
            intersection = Vector3.Lerp(closestLine2, closestLine2, 0.5f);
        }

        return intersection;
    }

I needed to calculate another plane perpendicular to an edge and get its line of intersection with another edge-plane. This line is intersected with the line calculated in the original post (instead of intersecting with a plane).

NOTE: Floating point inaccuracy means that the two lines don't actually intersect so I get the closest point on both lines and calculate a point halfway between them

Input:

a = (0.0, 3.4, 7.3), b = (-1.0, -1.7, -1.8), c = (1.0, -4.5, -1.0), d = (0.0, 1.4, -10.6)

Calculated Values:

planeABNormal = (-1.0, -5.1, -9.1), planeABPoint = (-0.5, 0.9, 2.8)

planeACNormal = (1.0, -7.9, -8.3), planeACPoint = (0.5, -0.6, 3.2)

planeBDNormal = (1.0, 3.1, -8.8), planeBDPoint = (-0.5, -0.2, -6.2)

planeCDNormal = (-1.0, 5.9, -9.6), planeCDPoint = (0.5, -1.6, -5.8)

linePoint1 = (1.0, -1.0, 3.6), lineDirection1 = (-29.6, -17.4, 13.0)

linePoint2 = (2.2, -3.1, -6.9), lineDirection2 = (22.2, 18.4, 9.0)

closestLine1 = (14.0, 6.7, -2.1), closestLine2 = (14.0, 6.7, -2.1)

intersection = (14.0, 6.7, -2.1) <- centre of circumsphere