Fernando D'Andrea
Fernando D'Andrea

Reputation: 307

How to get an hypersphere center point out of 5 surface points?

I'm trying to implement a voronoi sampler in 4-dimensional space, which should be a quite simple task, but picking a good starting points distribution, which can sometimes generate some quite skewed voronoi cells due to the lack of uniformity in distribution (I accept suggestions on how to generate uniform random point scattering).

So I started studying the case of generating centroidal voronoi cells and stumbled in the problem on the title. I'm starting with the assumption one can define a 4D hypersphere out of five surface points.

I went as far as getting some code to produce a 3D sphere center point out of three points in 3D space (in this case, the points are always in the sphere's equator), which I shared here. I also managed to find code that fives a sphere's center for tridimensional space out of four surface points, but it's got insanely large with lots of matrices determinant calculations, to which I admit it kind of crushed my hopes of extending it to 4D myself.

I found lot's of demonstrations in https://math.stackexchange.com, but that's not something I could readily transformable into code.

UPDATE UPDATE UPDATE!!!

Hello all, I'm finally on my weekend so I can give these some shots.

First, I arrived at the following C# (Unity3D) code I can visually confirm to be working wonders. It obviously fails when all points are coplanar, but that's pretty much expected.

static public Vector3 TetrahedronCircumcenter3D(Vector3 a, Vector3 b, Vector3 c, Vector3 d)
{
    Vector3 ba = b - a; //all points from an 'a' point of view
    Vector3 ca = c - a;
    Vector3 da = d - a;
    Vector3 crosscd = Vector3.Cross(ca, da); //perpendicular vectors to those above
    Vector3 crossdb = Vector3.Cross(da, ba);
    Vector3 crossbc = Vector3.Cross(ba, ca);
    return a + (
        (
            ba.sqrMagnitude * crosscd +
            ca.sqrMagnitude * crossdb +
            da.sqrMagnitude * crossbc
        ) *
        (0.5f / (ba.x * crosscd.x + ba.y * crosscd.y + ba.z * crosscd.z)) // half point
    );
}

As I said, it visually works.

I got this from a link I unfortunately lost that explains how to calc the circumcenter of a simplex from the crosspoint of perpendicular lines. These two words are important, for that was what I was looking for all the time: "simplex circumcenter". A simplex is the simplest shape one can get in R^n: a triangle in R^2, a tetrahedron in R^3 and a pentatope in R^4 and so forth.

I then decided to try to extend it to simply extend it to 4D using good sense. I first stumbled in the problem one can't expect to have a perpendicular cross-product between two vectors in anything but 3-D and (weird) 7-D, and I thought I was done. But a remark gave one next clue (which is kind of obvious, after one says it: you can get a perpendicular 4-D vector from the cross-product between 3 vectors in 4-D. Searching for the formula, I found something even better: code!

It is trivial to extend that code to work with Vector4.

The extended code looks like this:

static public Vector4 PentatopeCircumcenter4D(Vector4 a, Vector4 b, Vector4 c, Vector4 d, Vector4 e)
{
    Vector4 ba = b - a; //all points from an 'a' point of view
    Vector4 ca = c - a;
    Vector4 da = d - a;
    Vector4 ea = e - a;
    Vector4 crosscde = CrossProduct(ca, da, ea); //perpendicular vectors
    Vector4 crossdeb = CrossProduct(da, ea, ba);
    Vector4 crossebc = CrossProduct(ea, ba, ca);
    Vector4 crossbcd = CrossProduct(ba, ca, da);
    return a +
        (
            (
                crosscde * ba.sqrMagnitude +
                crossdeb * ca.sqrMagnitude +
                crossebc * da.sqrMagnitude +
                crossbcd * ea.sqrMagnitude
            ) *
            (0.5f / (ba.x * crosscde.x + ba.y * crosscde.y + ba.z * crosscde.z + ba.w * crosscde.w))
        );
}

Which begs the next question: how to confirm this is working?

Edit: the rationale is explained here: https://ctools.ece.utah.edu/Triangulation/TriangulationSphereCntr.pdf

Edit: It's not working in 4D. The distance to the calculated center from points a, b, c, d, and e is not the same. It's concordant for a, b, d. Points c and e show a different distance.

Edit: not having a linear algebra solver at hand in runtime and not being able to reduce the systems as appointed in other sources, I had to go the way of matrices. The code got a bit extensive, but nothing terrible. There's some 5x5 Matrix determinants in the way, but I got it working by checking distance between all points and the given center and it seems perfect now. The above tentative solutions do not work, though.

Upvotes: 0

Views: 267

Answers (1)

MBo
MBo

Reputation: 80197

Equation of hypersphere through determinant (like (29) here but with 4 dimensions) is quite gigantic for 4D.

There is a method that requires calculation of inverse matrix (rather simple process, standard problem in matrix algebra).

Choose one point as base X0 and make 4 vectors V1 = X1 - X0, V2 = X2 - X0 and so on. Calculate their half-lengths

L1 = V1.norm / 2 ...

Get normalized (unit) vectors:

U1 = V1 / (2*L1) ...

Solve equation to get radius-vector R (from X0 to hypersphere center)

 | R.x |     |  U1.x  U1.y  U1.z  U1.w |               | L1 |
 | R.y |     |  U2.x  U2.y  U2.z  U2.w |.Inversed *    | L2 |
 | R.z |  =  |  U3.x  U3.y  U3.z  U3.w |               | L3 |
 | R.w |     |  U4.x  U4.y  U4.z  U4.w |               | L4 |

And hypersphere center is

 C = X0 + R

Nice rationale for this equation is described here

It uses projection of radius-vector onto V-vectors, just like method of middle perpendiculars for 2D case of circumscribed circle around triangle

enter image description here enter image description here

Upvotes: 4

Related Questions