dennis2society.de » Blog Archive » Painless Tetrahedral Barycentric Mapping in C++

Painless Tetrahedral Barycentric Mapping in C++

2011/08/14 – 18:31

This article provides a set of functions to generate barycentric coordinates for surface mesh vertices in relation to a tetrahedral mesh. Additionally I provide some more functions to find the closest tetrahedron to a certain point P. As barycentric coordinates are invariant to translation, rotation and scaling, barycentric mapping is especially useful for efficiently updating a surface mesh according to an underlying deformable tetrahedral mesh.

This article was born out of frustration as I could not find any comprehensive guide on how to do this and I didn’t just want to “borrow” someone else’s code without really understanding what it actually does ( e.g. from SOFA Framework or PhysX SDK). It is primarily meant as an extended reminder for myself, but might be useful for other people dealing with tetrahedral meshes. This article will mostly contain a lot of code and only very few theoretical background and no mathematical proofs! At least that’s what I regard as “painless”. The advantage of my approach is that it simply relies on Vec3f and Vec4f types and avoids funny constructions like “MeshHash” as seen in the PhysX SDK. This will allow easy adaption to different VecXf type implementations and simplify integration into different applications.

For best understanding you should at least be familiar with the following concepts:

C / C++
surface meshes
tetrahedral meshes including useful data structures
basic knowledge of 3D vector math

Barycentric coordinate for a single vertex and a single tetrahedron

Calculating the barycentric coordinate for a single vertex in relation to a tetrahedron is actually pretty easy and involves calculating some determinants for 4×4 matrices:

/**
 * A Tetrahedron consisting of four vertices.
 */
struct Tetrahedron
{
    Vec3f v[4];
};
 
/**
 * Calculate the determinant for a 4x4 matrix based on this example:
 * http://www.euclideanspace.com/maths/algebra/matrix/functions/determinant/fourD/index.htm
 * This function takes four Vec4f as row vectors and calculates the resulting matrix' determinant
 * using the Laplace expansion.
 *
 */
const float Determinant4x4( const Vec4f& v0,
                            const Vec4f& v1,
                            const Vec4f& v2,
                            const Vec4f& v3 )
{
float det = v0.w*v1.z*v2.y*v3.x - v0.z*v1.w*v2.y*v3.x -
            v0.w*v1.y*v2.z*v3.x + v0.y*v1.w*v2.z*v3.x +
 
            v0.z*v1.y*v2.w*v3.x - v0.y*v1.z*v2.w*v3.x -
            v0.w*v1.z*v2.x*v3.y + v0.z*v1.w*v2.x*v3.y +
 
            v0.w*v1.x*v2.z*v3.y - v0.x*v1.w*v2.z*v3.y -
            v0.z*v1.x*v2.w*v3.y + v0.x*v1.z*v2.w*v3.y +
 
            v0.w*v1.y*v2.x*v3.z - v0.y*v1.w*v2.x*v3.z -
            v0.w*v1.x*v2.y*v3.z + v0.x*v1.w*v2.y*v3.z +
 
            v0.y*v1.x*v2.w*v3.z - v0.x*v1.y*v2.w*v3.z -
            v0.z*v1.y*v2.x*v3.w + v0.y*v1.z*v2.x*v3.w +
 
            v0.z*v1.x*v2.y*v3.w - v0.x*v1.z*v2.y*v3.w -
            v0.y*v1.x*v2.z*v3.w + v0.x*v1.y*v2.z*v3.w;
    return det;
}
 
/**
 * Calculate the actual barycentric coordinate from a point p0_ and the four 
 * vertices v0_ .. v3_ from a tetrahedron.
 */
const Vec4f GetBarycentricCoordinate( const Vec3f& v0_,
                                      const Vec3f& v1_,
                                      const Vec3f& v2_,
                                      const Vec3f& v3_,
                                      const Vec3f& p0_)
{
    Vec4f v0(v0_, 1);
    Vec4f v1(v1_, 1);
    Vec4f v2(v2_, 1);
    Vec4f v3(v3_, 1);
    Vec4f p0(p0_, 1);
    Vec4f barycentricCoord = Vec4f();
    const float det0 = Determinant4x4(v0, v1, v2, v3);
    const float det1 = Determinant4x4(p0, v1, v2, v3);
    const float det2 = Determinant4x4(v0, p0, v2, v3);
    const float det3 = Determinant4x4(v0, v1, p0, v3);
    const float det4 = Determinant4x4(v0, v1, v2, p0);
    barycentricCoord.x = (det1/det0);
    barycentricCoord.y = (det2/det0);
    barycentricCoord.z = (det3/det0);
    barycentricCoord.w = (det4/det0);
    return barycentricCoord;
}
 
const Vec4f GetBarycentricCoordinate( const Tetrahedron& t_,
                                      const Vec3f& p0_)
{
    return GetBarycentricCoordinate(t_.v[0], t_.v[1], t_.v[2], t_.v[3], p0_);
}

Finding the best tetrahedron for a vertex

In a real application we will not only work with a single vertex and only one tetrahedron. Instead there will be a list of tetrahedra and a list of surface-vertices and we will have to find the closest tetrahedron for each vertex. For each surface vertex P we have to iterate over all tetrahedra and either find the tetrahedron that contains our vertex or (if the vertex is not inside the tetrahedral mesh) find the tetrahedron whose centroid has the smallest distance to our vertex.

Check if the vertex P is inside the tetrahedron T. One easy way to do this once again involves calculating the determinants for five 4×4 matrices which can be constructed easily from our point P and the four vertices in the tetrahedron T.

/* Checks whether the specified point is inside the tetrahedron (by index)
*  based on this approach: http://steve.hollasch.net/cgindex/geometry/ptintet.html
*  @return true if the point_ is inside the tetrahedron (or on one of the four triangles), otherwise false
*  @param point_ the Vec3f point to be checked
*  @param tetra_ the tetrahedron
*/
const bool CheckPointInTetra(const Vec3f& point_, const Tetrahedron& tetra_)
{
    Vec4f v0(tetra_.v[0], 1);
    Vec4f v1(tetra_.v[1], 1);
    Vec4f v2(tetra_.v[2], 1);
    Vec4f v3(tetra_.v[3], 1);
    Vec4f p0(point_, 1);
    const float det0 = Determinant4x4(v0, v1, v2, v3);
    const float det1 = Determinant4x4(p0, v1, v2, v3);
    const float det2 = Determinant4x4(v0, p0, v2, v3);
    const float det3 = Determinant4x4(v0, v1, p0, v3);
    const float det4 = Determinant4x4(v0, v1, v2, p0);
    /**
    If by chance the Determinant det0 is 0, then your tetrahedron is degenerate (the points are coplanar).
    If any other Di=0, then P lies on boundary i (boundary i being that boundary formed by the three points other than Vi).
    If the sign of any Di differs from that of D0 then P is outside boundary i.
    If the sign of any Di equals that of D0 then P is inside boundary i.
    If P is inside all 4 boundaries, then it is inside the tetrahedron.
    As a check, it must be that D0 = D1+D2+D3+D4.
    The pattern here should be clear; the computations can be extended to simplicies of any dimension. (The 2D and 3D case are the triangle and the tetrahedron).
    If it is meaningful to you, the quantities bi = Di/D0 are the usual barycentric coordinates.
    Comparing signs of Di and D0 is only a check that P and Vi are on the same side of boundary i.
    */
    if (det0 != 0)
    {
        if (det0 < 0)
        {
            if ((det1 < 0) && (det2 < 0) && (det3 < 0) && (det4 < 0))
            {
                return true;
            }
        }
        if (det0 > 0)
        {
            if ((det1 > 0) && (det2 > 0) && (det3 > 0) && (det4 > 0))
            {
                return true;
            }
        }
    }
    return false;
}

If the point is not inside the tetrahedron we have to find the tetrahedron whose centroid has the smallest distance to our point P. The following code demonstrates how to calculate the distance between a point and a tetrahedron’s centroid.

/**
 * A Tetrahedron consisting of four vertices.
 */
struct Tetrahedron
{
    Vec3f v[4];
};
 
/**
 * Returns the centroid for the tetrahedron consisting of the four vertices
 */
const Vec3f GetTetrahedronCentroid( const Vec3f& v0_,
                                    const Vec3f& v1_,
                                    const Vec3f& v2_,
                                    const Vec3f& v3_)
{
    return (v0_ + v1_ + v2_ +v3_) / 4;
}
 
/**
 * Returns the centroid for the tetrahedron consisting of the four vertices.
 */
const Vec3f GetTetrahedronCentroid( const Tetrahedron& tet)
{
    Vec3f c;
    for (unsigned int i=0; i<4; ++i)
    {
        c += tet.v[i];
    }
    return c / 4;
}
 
/**
 * Returns the squared distance between a point and a tetrahedron's centroid.
 * This is useful when we don't need to know the exact distance but instead want to
 * compare between different distances, e.g. for finding the closest tetrahedron
 * to our point p. This also avoids calculating sqrt.
 * @return   the squared distance between point p and centroid of t
 * @param    p   the point
 * @param    t   the tetrahedron
 */
const float GetSquaredDistancePointToTetrahedron(const Vec3f& p, const Tetrahedron& t)
{
    const Vec3f centroid = GetTetrahedronCentroid(t);
    return (centroid - p).squaredLength();
}
 
/**
 * Returns the true distance between a point and a tetrahedron's centroid.
 */
const float GetDistancePointToTetrahedron(const Vec3f& p, const Tetrahedron& t)
{
    return sqrt(GetSquaredDistancePointToTetrahedron(p, t));
}

Algorithm

Given two lists (e.g. std::vector) for vertices and tetrahedra, the above functions can finally be put together into two simple functions. This combination also avoids iterating over all tetrahedra twice when a vertex is not inside a tetrahedron. This is important when your surface mesh is larger (in terms of volume) than the underlying tetrahdral mesh. This is done by calculating the distance between vertex P and tetrahedron T in the same step and checking whether the current distance is smaller than the previous one. Finally, if there is no tetrahedron containing the vertex, we will use the barycentric coord from the closest tetrahedron. Note that there is another list tetraIndices for tetrahdron-indices to store to which tetrahedron the vertex will be mapped.

 
std::vector<Vec4f> barycentricCoords; // the list of barycentric coordinates
std::vector<unsigned int> tetraIndices; // the list of tetrahedron indices
 
/**
 * Find the barycentric coordinate for a single vertex in a list of tetrahedra.
 */
void FindBarycentricCoordinateForVertex(const Vec3f& p, std::vector<Tetrahedron> tetras)
{
    float minDist = FLT_MAX;
    unsigned int closestTetrahedronIndex;
    for (unsigned int i = 0; i<tetras.size(); ++i)
    {
        if (CheckPointInTetra(p, tetras[i]))
        {
            tetraIndices.push_back(i);
            barycentricCoords.push_back(GetBarycentricCoordinate(tetras[i], p));
            return;
        }
        else
        {
            float dist = GetDistancePointToTetrahedron(p, tetras[i]);
            if (dist < minDist)
            {
                minDist = dist;
                closestTetrahedronIndex = i;
            }
        }
    }
    tetraIndices.push_back(closestTetrahedronIndex);
    barycentricCoords.push_back(GetBarycentricCoordinate(tetras[closestTetrahedronIndex], p));
}
 
/**
 * Find the tetrahedral barycentric coordinates for the surface vertices list.
 */
void FindBarycentricCoordinates(std::vector<Vec3f> surfaceVertices, std::vector<Tetrahedron> tetras)
{
    std::vector<Vec3f>::const_iterator vi;
    for (vi = surfaceVertices.begin(); vi<surfaceVertices.end(); ++vi)
    {
        FindBarycentricCoordinateForVertex(*vi, tetras);
    }
}

Some remarks

The Tetrahedron struct is for demonstration purposes only! In a real application the Tetrahedron struct will consist of four unsigned int values containing the indices in a tetrahedra-vertices array. This is helpful as several tetrahedra will share certain vertices. For this reason there exist some supplementary functions that accept single vertices (Vec3f) instead of the custom Tetrahedron struct. The functions for finding the barycentric coordinates in the first place have to be modified by yourself, though.

Applying the barycentric coordinates to update the surface mesh according to the tetrahedral mesh can be accomplished by simply doing the following. This will return the new position for the surface point mapped by the barycentric coordinate when the tetrahadron has moved, rotated, scaled or has been deformed (Screenshot on the right side).

void UpdateVertices( std::vector<Vec3f>& surfaceVertices, 
                     const std::vector<Tetrahedron>& tetras, 
                     const std::vector<unsigned int>& tetraIndices, 
                     const std::vector<Vec4f>& barycentricCoords)
{
    for (unsigned int vi = 0; vi<surfaceVertices.size(); ++vi)
    {
        surfaceVertices[vi] = GetPositionByBarycentricCoord(tetras[tetraIndices[vi]], barycentricCoords[vi], tetraIndices[vi]);
    }
    // generate surface normals
    // draw
    // etc.
}
 
const Vec3f GetPositionByBarycentricCoord( const Tetrahedron& t, 
                                           const Vec4f& barycentricCoord, 
                                           const unsigned int tetraIndex)
{
    Vec3f pos;
    for (unsigned int i=0; i<4; ++i)
    {
        pos += t.v[i] * barycentricCoord[i];
    }
    return pos;
}