Sphere::FromTetrahedron

  Sphere Sphere::FromTetrahedron(const glm::vec3& p1, const glm::vec3& p2,
                                 const glm::vec3& p3, const glm::vec3& p4)
  {
    assert(AreEqual(p1, p2) == false);
    assert(AreEqual(p1, p3) == false);
    assert(AreEqual(p1, p4) == false);
    assert(AreEqual(p2, p3) == false);
    assert(AreEqual(p2, p4) == false);
    assert(AreEqual(p3, p4) == false);

    glm::mat4 matrix(1.0f);
    matrix = glm::row(matrix, 0, glm::vec4(p1.x, p2.x, p3.x, p4.x));
    matrix = glm::row(matrix, 1, glm::vec4(p1.y, p2.y, p3.y, p4.y));
    matrix = glm::row(matrix, 2, glm::vec4(p1.z, p2.z, p3.z, p4.z));
    matrix = glm::row(matrix, 3, glm::vec4(1.0f, 1.0f, 1.0f, 1.0f));

    // Function that assert if determinant is 0
    const float a = Determinant(matrix);

    glm::mat4 D = matrix;
    glm::vec3 center(0.0f);

    const glm::vec4 squares = {
      GetSquaredLength(p1),
      GetSquaredLength(p2),
      GetSquaredLength(p3),
      GetSquaredLength(p4),
    };

    // D[0] = squares;
    D        = glm::row(D, 0, squares);
    center.x = glm::determinant(D);

    // D[1] = matrix[0];
    D        = glm::row(D, 1, glm::row(matrix, 0));
    center.y = -glm::determinant(D);

    // D[2] = matrix[1];
    D        = glm::row(D, 2, glm::row(matrix, 1));
    center.z = glm::determinant(D);

    center /= (2.0f * a);
    return Sphere(center, glm::distance(p1, center));