9
\$\begingroup\$

I'm looking for the most minimized equation to find the center coordinates and the radius of a tetrahedron circumsphere given four 3D points.

What I found on the internet mainly deal with the circum sphere of a flat 3D triangle, or some rough mathematical definitions, or some very single case such as regular tetrahedrons. Anyway I managed to find the equation below but I missed something :

    ->  ->      ->
let d1, d2, and d3 three vectors of any face of the triangle :

    | d1x  d1y  d1z |   | x |   | d1^2 |
2 * | d2x  d2y  d2z | * | y | = | d2^2 |
    | d3x  d3y  d3z |   | z |   | d3^2 |

My knowledge in this field has its limits but I think I can handle matrices and vector operations. But is the right part of the equation the square of the norm of each vectors ? (which are into a vector). Is the equation valid ? Is it just the writer who lazely forgot to write |d1|^2 ? Or Is it a common way to define some mathematical property.

PS : It's for a Delaunay Triangulation implementation. The equation (number 9) is in the following link : https://www2.mps.mpg.de/homes/daly/CSDS/t4h/tetra.htm

\$\endgroup\$
  • 4
    \$\begingroup\$ Try maths stackexchange. \$\endgroup\$ – Majte Oct 24 '15 at 2:10
  • \$\begingroup\$ Thanks I've found a way to compute the circumsphere there ! \$\endgroup\$ – herme5 Oct 24 '15 at 3:07
  • 1
    \$\begingroup\$ @JamesAMD the link is www2.mps.mpg.de/homes/daly/CSDS/t4h/tetra.htm. \$\endgroup\$ – herme5 Oct 24 '15 at 11:02
  • 3
    \$\begingroup\$ @herme5, feel free to post your own answer here on how you are calculating the answer. Many people may come here in the future hoping to find the answer, and you sharing it will be valuable to them. It's totally acceptable to post your own answer and even accept it. \$\endgroup\$ – Tim Holt Feb 4 '18 at 23:57
  • 2
    \$\begingroup\$ Thanks for the notice @TimHolt. I’ll do it ! Nevertheless I am not sure anymore how I did it, it was more than 2 years ago ! just let me find and take a look at my old implementation \$\endgroup\$ – herme5 Feb 5 '18 at 1:53
2
\$\begingroup\$

While this is an ancient thread, I thought it may be nice for posterity to have a bit to reference. The source of the formula is from Geometric Tools for Computer Graphics by Philip J. Schneider and David H. Eberly. Something to note, according to the text

The tetrahedron V0, V1, V2, V3 is ordered so that it is isomorphic to the canonical one (0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1) .

As I understand isomorphism, there can be several different meanings when used in geometry. If he means isomorphic with regard to graph theory, then the following code should behave correctly, as the topology of any tetrahedron is the same (K4, a complete graph). I tested the results of the function against wolfram alpha using various permutations in the ordering of the canonical vertices, and I saw no difference in the outcome. If the ordering proves to be a problem, I suggest examining the normal of the triangle formed by the vertices V1,V2,V3 upon input to this function, and treating the points like a half-space with a dot-product test to figure out if that triangle is facing the right way. If it's not, a simple std::swap of any two of the triangle's vertices will reverse the direction of the normal and you may continue. But like I said, I saw no difference with various permutations.

Here is the translated code without using matrices to avoid any implementation confusion, it's fairly straight-forward;

void Circumsphere(const Vec3& v0, const Vec3& v1, const Vec3& v2, const Vec3& v3, Vec3* center, float* radius)
{
  //Create the rows of our "unrolled" 3x3 matrix
  Vec3 Row1 = v1 - v0;
  float sqLength1 = length2(Row1);
  Vec3 Row2 = v2 - v0;
  float sqLength2 = length2(Row2);
  Vec3 Row3 = v3 - v0;
  float sqLength3 = length2(Row3);

  //Compute the determinant of said matrix
  const float determinant =   Row1.x * (Row2.y * Row3.z - Row3.y * Row2.z)
                            - Row2.x * (Row1.y * Row3.z - Row3.y * Row1.z)
                            + Row3.x * (Row1.y * Row2.z - Row2.y * Row1.z);

  // Compute the volume of the tetrahedron, and precompute a scalar quantity for re-use in the formula
  const float volume = determinant / 6.f;
  const float iTwelveVolume = 1.f / (volume * 12.f);

  center->x = v0.x + iTwelveVolume * ( ( Row2.y * Row3.z - Row3.y * Row2.z) * sqLength1 - (Row1.y * Row3.z - Row3.y * Row1.z) * sqLength2 + (Row1.y * Row2.z - Row2.y * Row1.z) * sqLength3 );
  center->y = v0.y + iTwelveVolume * (-( Row2.x * Row3.z - Row3.x * Row2.z) * sqLength1 + (Row1.x * Row3.z - Row3.x * Row1.z) * sqLength2 - (Row1.x * Row2.z - Row2.x * Row1.z) * sqLength3 );
  center->z = v0.z + iTwelveVolume * ( ( Row2.x * Row3.y - Row3.x * Row2.y) * sqLength1 - (Row1.x * Row3.y - Row3.x * Row1.y) * sqLength2 + (Row1.x * Row2.y - Row2.x * Row1.y) * sqLength3 );

  //Once we know the center, the radius is clearly the distance to any vertex
  *radius = length(*center - v0);
}
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.