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
May have figured it out (tests look good so far):
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