Barycentric coordinates of a tetrahedron

Question

I would like to ask for help regarding barycentric coordinates of a tetrahedron:

Following an approach I found here: http://www.cdsimpson.net/2014/10/barycentric-coordinates.html i implemented a C++ function for finding barycentric coordinates of a point in a tetrahedron:

float ScTP(const Vec &a, const Vec &b, const Vec &c)
{
    // computes scalar triple product
    return Dot(a, Cross(b, c));
}

Vec4f bary_tet(Vec3f a, Vec3f b, Vec3f c, Vec3f d, Vec3f p)
{
    float va, vb, vc, vd, v;
    Vec3f vap = p - a;
    Vec3f vbp = p - b;
    Vec3f vcp = p - c;
    Vec3f vdp = p - d;

    Vec3f vab = b - a;
    Vec3f vac = c - a;
    Vec3f vad = d - a;

    Vec3f vbc = c - b;
    Vec3f vbd = d - b;
    // ScTP computes the scalar triple product
    va = ScTP(vbp, vbd, vbc) * 1 / 6;
    vb = ScTP(vap, vac, vad) * 1 / 6;
    vc = ScTP(vap, vad, vab) * 1 / 6;
    vd = ScTP(vap, vab, vac) * 1 / 6;
    v = 1 / ScTP(vab, vac, vad) * 1 / 6;
    return Vec4f(va*v, vb*v, vc*v, vd*v);
}

However, my code seems to calculate slightly wrong barycentric coordinates - comparing my results with a reference implementation from here: http://dennis2society.de/painless-tetrahedral-barycentric-mapping my four barycentric values are each smaller the values calculated by the reference implementation.

Does anyone spot any error in my implementation? Many thanks for help!

Snapper

how much smaller? Couldn't it be just float vs double precision? (didn't run your code, as it looks not very complete, ScTP looks missing) (also you shouldn't copy input parameters like this... use const Vec3f & a rather, if the performance is important (probably is, because otherwise why would you use float instead of double) – Grate

unfortunately not, my values are roughly 1/40 smaller.. – Snapper

thanks for your suggestion, i did not yet look for speed, but i will add it for the final implementation. also I added the ScTP function to the code – Snapper

I'm not sure, but I think "v" may be written other way you intended. You want (1/6)*ScTP and that to flip? But you do 1/ScTP and then /6 further. Btw, you can simply drop all the 1/6 if I didn't lost all my math skills completely, if you use only va/v, ... – Grate

Answer 1

Blind guess:

Vec4f bary_tet(const Vec3f & a, const Vec3f & b, const Vec3f & c, const Vec3f & d, const Vec3f & p)
{
    Vec3f vap = p - a;
    Vec3f vbp = p - b;

    Vec3f vab = b - a;
    Vec3f vac = c - a;
    Vec3f vad = d - a;

    Vec3f vbc = c - b;
    Vec3f vbd = d - b;
    // ScTP computes the scalar triple product
    float va6 = ScTP(vbp, vbd, vbc);
    float vb6 = ScTP(vap, vac, vad);
    float vc6 = ScTP(vap, vad, vab);
    float vd6 = ScTP(vap, vab, vac);
    float v6 = 1 / ScTP(vab, vac, vad);
    return Vec4f(va6*v6, vb6*v6, vc6*v6, vd6*v6);
}

Grate

many thanks, your version works!! and i understand the error in my code. – Snapper