5

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!

christianlehmann
  • 95
  • 2
  • 9
  • 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`) – Ped7g Jul 23 '16 at 19:39
  • unfortunately not, my values are roughly 1/40 smaller.. – christianlehmann Jul 23 '16 at 19:40
  • 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 – christianlehmann Jul 23 '16 at 19:47
  • 1
    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, ... – Ped7g Jul 23 '16 at 20:16

1 Answers1

9

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);
}
Ped7g
  • 16,236
  • 3
  • 26
  • 63