я уверен, что вы все затаили дыхание и ждёте, когда же я вывалю кучу кода, которая покажет, как эта вся херня реально работает. без проблем.
итак. самую важную часть рэйкастера я уже описал . остальное есть в статье, но фиг с ним, приведу полный код:
public Raycast!VT gjkraycast(bool checkRayStart=true, int maxiters=32, double distEps=0.0001, VT, CT) (auto ref CT abody, in auto ref VT rayO, in auto ref VT rayD) if (IsGoodVertexStorage!(CT, VT)) {
Raycast!VT res;
immutable VT start = rayO;
VT r = rayD;
VT.Float maxlen = r.length;
r.normalize;
VT n; // normal at the hit point
VT x = start; // current closest point on the ray
VT.Float lambda = 0;
VT a = VT.Invalid, b = VT.Invalid; // simplex
VT v = x-abody.centroid;
VT.Float distsq = VT.Float.infinity;
int itersLeft = maxiters;
while (itersLeft > 0) {
VT p = abody.support(v);
VT w = x-p;
VT.Float dvw = v.dot(w);
if (dvw > 0) {
VT.Float dvr = v.dot(r);
if (dvr >= -(EPSILON!(VT.Float)*EPSILON!(VT.Float))) return res;
lambda = lambda-dvw/dvr;
if (lambda > maxlen) return res;
x = start+r*lambda;
n = v;
}
// reduce simplex
if (a.valid) {
if (b.valid) {
VT p1 = x.projectToSeg(a, p);
VT p2 = x.projectToSeg(p, b);
if (p1.distanceSquared(x) < p2.distanceSquared(x)) {
b = p;
distsq = p1.distanceSquared(x);
} else {
a = p;
distsq = p2.distanceSquared(x);
}
VT ab = b-a;
VT ax = x-a;
v = ab.tripleProduct(ax, ab);
} else {
b = p;
VT ab = b-a;
VT ax = x-a;
v = ab.tripleProduct(ax, ab);
}
} else {
a = p;
v = -v;
}
if (distsq <= cast(VT.Float)distEps) break;
--itersLeft;
}
if (itersLeft < 1) return res; // alas, out of iterations
// result
res.p = x;
res.n = n.normalized;
res.dist = lambda;
return res;
}
да, это всё. практически один-в-один код из статьи.
теперь перейдём к sweep tests. повторю: самый простой метод реализовать это — «надуть» один полигон на размер другого; тогда «другой» полигон можно редуцировать в точку, и сделать обычный рэйкаст. так работает AABB vs AABB sweep test, так работает проверка столкновений в кваке (и думе, по-моему; извините, точно не помню). методов «надувания» полигонов есть разных, но в случае GJK нам интересен один: minkowski difference.
тут есть нюанс: если minkowski sum пишут везде одинаково: A+B, то с minkowski difference есть непонятки. то это A-B, то это -A+B, то -A-B, то ещё какая-то ебанина. на самом деле это всё одно и то же, просто операция вычитания некоммутативна (очевидно), и её пишут в зависимости от того, как у них там что представлено. в нашем случае это будет вот так:
VT p = bbody.support(v)-(abody.support(-v)-abody.centroid);
да-да. это все важные изменения в рэйкастере! узрите! здесь у нас abody летит с линейной скоростью lvelA, и мы смотрим, долбанётся ли оно об bbody.
public Raycast!VT gjksweep(bool checkRayStart=true, int maxiters=32, double distEps=0.0001, VT, CT) (auto ref CT abody, auto ref CT bbody, in auto ref VT lvelA) {
Raycast!VT res;
immutable VT start = abody.centroid; //(!!!) trace from abody
VT r = lvelA; //(!!!) relative motion
VT.Float maxlen = r.length;
r.normalize;
VT n; // normal at the hit point
VT x = start; // current closest point on the ray
VT.Float lambda = 0;
VT a = VT.Invalid, b = VT.Invalid; // simplex
VT v = x-bbody.centroid; //(!!!)
VT.Float distsq = VT.Float.infinity;
int itersLeft = maxiters;
while (itersLeft > 0) {
VT p = bbody.support(v)-(abody.support(-v)-abody.centroid); //(!!!)
VT w = x-p;
VT.Float dvw = v.dot(w);
if (dvw > 0) {
VT.Float dvr = v.dot(r);
if (dvr >= -(EPSILON!(VT.Float)*EPSILON!(VT.Float))) return res; //(***)
lambda = lambda-dvw/dvr;
if (lambda > maxlen) return res;
x = start+r*lambda;
n = v;
}
// reduce simplex (all the following code is not changed)
все изменённые строки помечены (!!!), и их, как видите, просто дофига.
самая важная изменённая строка — последняя. именно в ней мы «надуваем» bbody. замечу, что для простоты все вершины у меня в world space: отсюда взялось -abody.centroid. если не привести координаты abody в object space, то bbody раздует так, что мало не покажется. нам же не нужно раздувать его на расстояние от (0,0) до abody!
ещё обратите внимание на строку с (***): это, на самом деле, всего лишь dvr >= 0. но если не впилить туда эпсилон, то из-за неточностей плавающей точки алгоритм в некоторых случаях прокручивает все итерации, прежде чем сдаться (когда луч проходит очень рядом с вершиной, и ещё при некоторых условиях).
могут спросить, почему лимит в 32 итерации. отвечаю: а просто так. на деле их обычно две-четыре, иногда чуть больше. и про distEps: ответ тот же самый. заметьте только, что distEps не должен быть слишком близок к эпсилону, а то итераций станет больше. реально эти числа взяты с потолка и впилены без особой теоретической базы — потому что работают. инженерная специфика.
