i > u

Painless Tetrahedral Barycentric Mapping in C++

2011-08-14 – 18:31

Screenshot of a liver surface mesh mapped to a tetrahedral mesh

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:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73

/**
 * 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.

    1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48

    /* 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.

    1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54

    /**
 * 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.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44

 
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);
    }
}

Barycentrically mapped torus surface meshSome 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).

    1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

    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;
}
References
http://www.euclideanspace.com/maths/algebra/matrix/functions/determinant/fourD/index.htm
BARYCENTRIC COORDINATES IN TETRAHEDRONS (PDF) – John Burkardt Information Technology Department Virginia Tech
http://steve.hollasch.net/cgindex/geometry/ptintet.html
Sofa Framework
PhysX SDK

  1. 5 Responses to “Painless Tetrahedral Barycentric Mapping in C++”

  2. Little confusion about mapping vertex to barycentric coordinates while updating .. since we are calculating barycentric coordinates with respect to the surface vertex, GetPositionByBarycentricCoord always giving the same position and vertex is not getting updated .. any help would be greatly appreciated..

    By jyothsna on Jun 7, 2013

  3. If I understand you correctly you have calculated the barycentric coordinates for surface mesh vertices in relation to the closest tetrahedron respectively in relation to the tetrahedron that contains the current surface vertex.

    Motion to the barycentrically mapped vertex will only happen when there occurs some rigid or non-rigid transformation to the parent tetrahedron (i.e. to the tetrahedron that contains your surface vertex), such as translation, rotation or deformation.
    To do this you have to multiply the Tetrahedron’s updated vertices (i.e. after the transformation) with the previously calculated barycentric coordinates for the current surface vertex (as described in GetPositionByBarycentricCoord). This will give you the new position of the surface vertex with respect to the transformed tetrahedron.
    Best.
    Dennis

    By dennis2society on Jun 7, 2013

  4. Dear Dennis

    Thank you very much for clarifying my doubt, with your explanation I was able to figure out that I was updating tetra info upfront before calculating barycenters. Now its all working fine, thank you very much for presenting this nice article..

    :) cheers..

    By Jyothsna on Jun 8, 2013

  5. Hi Dennis,

    Nice article!
    I have a similar problem to solve; I need to write a C++ function that will test if a point in 3D is inside a square pyramid (5 vertices). Can you please let me know how your code could be modified for this? Also do you think that this method is more efficient than trying to compute face normals?

    Cheers

    By Fabien on Aug 15, 2013

  6. Hi Fabien.
    One option for the point-inside-5-point-pyramid test is to simply sub-divide your pyramid into two tetrahedra and do the point-in-tetra-test twice (once for each resulting tetrahedron).
    Sub-dividing the pyramid into tetrahedra is not too complicated and mostly depends on the vertex-order of your pyramid.
    Another approach would be a point-in-frustum-test as seen here:
    http://www.lighthouse3d.com/tutorials/view-frustum-culling/
    This could easily be extended to a true pyramid instead of a frustum. This requires calculating the normals of the frustum planes and some dinstances.

    However, I have no idea which of these methods would be the more efficient…

    By dennis2society on Aug 16, 2013

Post a Comment

Insert a smiley:


To prevent spam and comments of people who don't
know about mathematics you have to enter the result
into the field on the right side.
Hint: Possible answers are '-3', '-2', '0', '1' and '6x'.
Some challenges have more than one correct answer!

Please note: Comments will be approved manually!
It may take a few hours to check if yours
is legitimate even after solving the captcha..
Spam protection

Captcha    

© 2002-2024