Chociaż jest to starożytny wątek, pomyślałem, że dla potomności może być coś ciekawego. Źródło formuły pochodzi z Geometric Tools for Computer Graphics autorstwa Philipa J. Schneidera i Davida H. Eberly. Coś do odnotowania, zgodnie z tekstem
Czworościan V0, V1, V2, V3 jest uporządkowany tak, że jest izomorficzny w stosunku do kanonicznego (0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1) ).
Jak rozumiem izomorfizm , w geometrii może istnieć kilka różnych znaczeń. Jeśli chodzi o izomorficzny w odniesieniu do teorii grafów, to poniższy kod powinien zachowywać się poprawnie, ponieważ topologia dowolnego czworościanu jest taka sama (K4, pełny wykres). I przetestowane wyniki funkcji przeciwko Wolfram Alpha stosując różne permutacje w zamawiania kanonicznych wierzchołków i widziałem żadnej różnicy w wynikach. Jeśli uporządkowanie okaże się problemem, sugeruję zbadanie normalności trójkąta utworzonego przez wierzchołki V1, V2, V3 po wejściu do tej funkcji i potraktowanie punktów jak półprzestrzeń testem iloczynu, aby ustalić jeśli ten trójkąt jest skierowany we właściwą stronę. Jeśli nie, to prostestd::swap
dowolnych dwóch wierzchołków trójkąta odwróci kierunek normalny i możesz kontynuować. Ale jak powiedziałem, nie widziałem różnicy przy różnych permutacjach.
Oto przetłumaczony kod bez użycia macierzy, aby uniknąć nieporozumień związanych z implementacją, jest dość prosty;
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);
}