Барицентрические координаты тетраэдра

Я хотел бы попросить помощи относительно барицентрических координат тетраэдра:

Следуя подходу, который я нашел здесь: http://www.cdsimpson.net/2014/10/barycentric-coordinates.html я реализовал функцию C++ для нахождения барицентрических координат точки в тетраэдре:

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

Однако мой код, кажется, вычисляет немного неправильные барицентрические координаты — сравнивая мои результаты с эталонной реализацией отсюда: http://dennis2society.de/painless-tetrahedral-barycentric-mapping каждое из моих четырех барицентрических значений меньше значений, рассчитанных эталонной реализацией.

Кто-нибудь заметил ошибку в моей реализации? Большое спасибо за помощь!


person christianlehmann    schedule 23.07.2016    source источник
comment
насколько меньше? Разве это не может быть просто float против double точности? (не запускал ваш код, так как он выглядит не очень полным, ScTP выглядит отсутствующим) (также вам не следует копировать входные параметры, подобные этому... скорее используйте const Vec3f & a, если важна производительность (вероятно, потому что иначе зачем бы вы используете float вместо double)   -  person Ped7g    schedule 23.07.2016
comment
к сожалению нет, мои значения примерно на 1/40 меньше..   -  person christianlehmann    schedule 23.07.2016
comment
спасибо за ваше предложение, я еще не искал скорость, но я добавлю ее для окончательной реализации. также я добавил в код функцию ScTP   -  person christianlehmann    schedule 23.07.2016
comment
Я не уверен, но я думаю, что v может быть написано иначе, чем вы предполагали. Хочешь (1/6)*ScTP и что флипнуть? Но вы делаете 1/ScTP и далее /6. Кстати, вы можете просто отбросить все 1/6, если я не потерял все свои математические навыки полностью, если вы используете только va/v, ...   -  person Ped7g    schedule 23.07.2016


Ответы (1)


Слепое предположение:

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);
}
person Ped7g    schedule 23.07.2016
comment
большое спасибо, ваша версия работает!! и я понимаю ошибку в моем коде. - person christianlehmann; 23.07.2016