diff options
Diffstat (limited to 'src/BulletCollision/NarrowPhaseCollision/btMprPenetration.h')
-rw-r--r-- | src/BulletCollision/NarrowPhaseCollision/btMprPenetration.h | 1768 |
1 files changed, 884 insertions, 884 deletions
diff --git a/src/BulletCollision/NarrowPhaseCollision/btMprPenetration.h b/src/BulletCollision/NarrowPhaseCollision/btMprPenetration.h index ea946e2da..358bc95d8 100644 --- a/src/BulletCollision/NarrowPhaseCollision/btMprPenetration.h +++ b/src/BulletCollision/NarrowPhaseCollision/btMprPenetration.h @@ -1,884 +1,884 @@ -
-/***
- * ---------------------------------
- * Copyright (c)2012 Daniel Fiser <danfis@danfis.cz>
- *
- * This file was ported from mpr.c file, part of libccd.
- * The Minkoski Portal Refinement implementation was ported
- * to OpenCL by Erwin Coumans for the Bullet 3 Physics library.
- * The original MPR idea and implementation is by Gary Snethen
- * in XenoCollide, see http://github.com/erwincoumans/xenocollide
- *
- * Distributed under the OSI-approved BSD License (the "License");
- * see <http://www.opensource.org/licenses/bsd-license.php>.
- * This software is distributed WITHOUT ANY WARRANTY; without even the
- * implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
- * See the License for more information.
- */
-
-///2014 Oct, Erwin Coumans, Use templates to avoid void* casts
-
-#ifndef BT_MPR_PENETRATION_H
-#define BT_MPR_PENETRATION_H
-
-#define BT_DEBUG_MPR1
-
-#include "LinearMath/btTransform.h"
-#include "LinearMath/btAlignedObjectArray.h"
-
-//#define MPR_AVERAGE_CONTACT_POSITIONS
-
-struct btMprCollisionDescription
-{
- btVector3 m_firstDir;
- int m_maxGjkIterations;
- btScalar m_maximumDistanceSquared;
- btScalar m_gjkRelError2;
-
- btMprCollisionDescription()
- : m_firstDir(0, 1, 0),
- m_maxGjkIterations(1000),
- m_maximumDistanceSquared(1e30f),
- m_gjkRelError2(1.0e-6)
- {
- }
- virtual ~btMprCollisionDescription()
- {
- }
-};
-
-struct btMprDistanceInfo
-{
- btVector3 m_pointOnA;
- btVector3 m_pointOnB;
- btVector3 m_normalBtoA;
- btScalar m_distance;
-};
-
-#ifdef __cplusplus
-#define BT_MPR_SQRT sqrtf
-#else
-#define BT_MPR_SQRT sqrt
-#endif
-#define BT_MPR_FMIN(x, y) ((x) < (y) ? (x) : (y))
-#define BT_MPR_FABS fabs
-
-#define BT_MPR_TOLERANCE 1E-6f
-#define BT_MPR_MAX_ITERATIONS 1000
-
-struct _btMprSupport_t
-{
- btVector3 v; //!< Support point in minkowski sum
- btVector3 v1; //!< Support point in obj1
- btVector3 v2; //!< Support point in obj2
-};
-typedef struct _btMprSupport_t btMprSupport_t;
-
-struct _btMprSimplex_t
-{
- btMprSupport_t ps[4];
- int last; //!< index of last added point
-};
-typedef struct _btMprSimplex_t btMprSimplex_t;
-
-inline btMprSupport_t *btMprSimplexPointW(btMprSimplex_t *s, int idx)
-{
- return &s->ps[idx];
-}
-
-inline void btMprSimplexSetSize(btMprSimplex_t *s, int size)
-{
- s->last = size - 1;
-}
-
-#ifdef DEBUG_MPR
-inline void btPrintPortalVertex(_btMprSimplex_t *portal, int index)
-{
- printf("portal[%d].v = %f,%f,%f, v1=%f,%f,%f, v2=%f,%f,%f\n", index, portal->ps[index].v.x(), portal->ps[index].v.y(), portal->ps[index].v.z(),
- portal->ps[index].v1.x(), portal->ps[index].v1.y(), portal->ps[index].v1.z(),
- portal->ps[index].v2.x(), portal->ps[index].v2.y(), portal->ps[index].v2.z());
-}
-#endif //DEBUG_MPR
-
-inline int btMprSimplexSize(const btMprSimplex_t *s)
-{
- return s->last + 1;
-}
-
-inline const btMprSupport_t *btMprSimplexPoint(const btMprSimplex_t *s, int idx)
-{
- // here is no check on boundaries
- return &s->ps[idx];
-}
-
-inline void btMprSupportCopy(btMprSupport_t *d, const btMprSupport_t *s)
-{
- *d = *s;
-}
-
-inline void btMprSimplexSet(btMprSimplex_t *s, size_t pos, const btMprSupport_t *a)
-{
- btMprSupportCopy(s->ps + pos, a);
-}
-
-inline void btMprSimplexSwap(btMprSimplex_t *s, size_t pos1, size_t pos2)
-{
- btMprSupport_t supp;
-
- btMprSupportCopy(&supp, &s->ps[pos1]);
- btMprSupportCopy(&s->ps[pos1], &s->ps[pos2]);
- btMprSupportCopy(&s->ps[pos2], &supp);
-}
-
-inline int btMprIsZero(float val)
-{
- return BT_MPR_FABS(val) < FLT_EPSILON;
-}
-
-inline int btMprEq(float _a, float _b)
-{
- float ab;
- float a, b;
-
- ab = BT_MPR_FABS(_a - _b);
- if (BT_MPR_FABS(ab) < FLT_EPSILON)
- return 1;
-
- a = BT_MPR_FABS(_a);
- b = BT_MPR_FABS(_b);
- if (b > a)
- {
- return ab < FLT_EPSILON * b;
- }
- else
- {
- return ab < FLT_EPSILON * a;
- }
-}
-
-inline int btMprVec3Eq(const btVector3 *a, const btVector3 *b)
-{
- return btMprEq((*a).x(), (*b).x()) && btMprEq((*a).y(), (*b).y()) && btMprEq((*a).z(), (*b).z());
-}
-
-template <typename btConvexTemplate>
-inline void btFindOrigin(const btConvexTemplate &a, const btConvexTemplate &b, const btMprCollisionDescription &colDesc, btMprSupport_t *center)
-{
- center->v1 = a.getObjectCenterInWorld();
- center->v2 = b.getObjectCenterInWorld();
- center->v = center->v1 - center->v2;
-}
-
-inline void btMprVec3Set(btVector3 *v, float x, float y, float z)
-{
- v->setValue(x, y, z);
-}
-
-inline void btMprVec3Add(btVector3 *v, const btVector3 *w)
-{
- *v += *w;
-}
-
-inline void btMprVec3Copy(btVector3 *v, const btVector3 *w)
-{
- *v = *w;
-}
-
-inline void btMprVec3Scale(btVector3 *d, float k)
-{
- *d *= k;
-}
-
-inline float btMprVec3Dot(const btVector3 *a, const btVector3 *b)
-{
- float dot;
-
- dot = btDot(*a, *b);
- return dot;
-}
-
-inline float btMprVec3Len2(const btVector3 *v)
-{
- return btMprVec3Dot(v, v);
-}
-
-inline void btMprVec3Normalize(btVector3 *d)
-{
- float k = 1.f / BT_MPR_SQRT(btMprVec3Len2(d));
- btMprVec3Scale(d, k);
-}
-
-inline void btMprVec3Cross(btVector3 *d, const btVector3 *a, const btVector3 *b)
-{
- *d = btCross(*a, *b);
-}
-
-inline void btMprVec3Sub2(btVector3 *d, const btVector3 *v, const btVector3 *w)
-{
- *d = *v - *w;
-}
-
-inline void btPortalDir(const btMprSimplex_t *portal, btVector3 *dir)
-{
- btVector3 v2v1, v3v1;
-
- btMprVec3Sub2(&v2v1, &btMprSimplexPoint(portal, 2)->v,
- &btMprSimplexPoint(portal, 1)->v);
- btMprVec3Sub2(&v3v1, &btMprSimplexPoint(portal, 3)->v,
- &btMprSimplexPoint(portal, 1)->v);
- btMprVec3Cross(dir, &v2v1, &v3v1);
- btMprVec3Normalize(dir);
-}
-
-inline int portalEncapsulesOrigin(const btMprSimplex_t *portal,
- const btVector3 *dir)
-{
- float dot;
- dot = btMprVec3Dot(dir, &btMprSimplexPoint(portal, 1)->v);
- return btMprIsZero(dot) || dot > 0.f;
-}
-
-inline int portalReachTolerance(const btMprSimplex_t *portal,
- const btMprSupport_t *v4,
- const btVector3 *dir)
-{
- float dv1, dv2, dv3, dv4;
- float dot1, dot2, dot3;
-
- // find the smallest dot product of dir and {v1-v4, v2-v4, v3-v4}
-
- dv1 = btMprVec3Dot(&btMprSimplexPoint(portal, 1)->v, dir);
- dv2 = btMprVec3Dot(&btMprSimplexPoint(portal, 2)->v, dir);
- dv3 = btMprVec3Dot(&btMprSimplexPoint(portal, 3)->v, dir);
- dv4 = btMprVec3Dot(&v4->v, dir);
-
- dot1 = dv4 - dv1;
- dot2 = dv4 - dv2;
- dot3 = dv4 - dv3;
-
- dot1 = BT_MPR_FMIN(dot1, dot2);
- dot1 = BT_MPR_FMIN(dot1, dot3);
-
- return btMprEq(dot1, BT_MPR_TOLERANCE) || dot1 < BT_MPR_TOLERANCE;
-}
-
-inline int portalCanEncapsuleOrigin(const btMprSimplex_t *portal,
- const btMprSupport_t *v4,
- const btVector3 *dir)
-{
- float dot;
- dot = btMprVec3Dot(&v4->v, dir);
- return btMprIsZero(dot) || dot > 0.f;
-}
-
-inline void btExpandPortal(btMprSimplex_t *portal,
- const btMprSupport_t *v4)
-{
- float dot;
- btVector3 v4v0;
-
- btMprVec3Cross(&v4v0, &v4->v, &btMprSimplexPoint(portal, 0)->v);
- dot = btMprVec3Dot(&btMprSimplexPoint(portal, 1)->v, &v4v0);
- if (dot > 0.f)
- {
- dot = btMprVec3Dot(&btMprSimplexPoint(portal, 2)->v, &v4v0);
- if (dot > 0.f)
- {
- btMprSimplexSet(portal, 1, v4);
- }
- else
- {
- btMprSimplexSet(portal, 3, v4);
- }
- }
- else
- {
- dot = btMprVec3Dot(&btMprSimplexPoint(portal, 3)->v, &v4v0);
- if (dot > 0.f)
- {
- btMprSimplexSet(portal, 2, v4);
- }
- else
- {
- btMprSimplexSet(portal, 1, v4);
- }
- }
-}
-template <typename btConvexTemplate>
-inline void btMprSupport(const btConvexTemplate &a, const btConvexTemplate &b,
- const btMprCollisionDescription &colDesc,
- const btVector3 &dir, btMprSupport_t *supp)
-{
- btVector3 seperatingAxisInA = dir * a.getWorldTransform().getBasis();
- btVector3 seperatingAxisInB = -dir * b.getWorldTransform().getBasis();
-
- btVector3 pInA = a.getLocalSupportWithMargin(seperatingAxisInA);
- btVector3 qInB = b.getLocalSupportWithMargin(seperatingAxisInB);
-
- supp->v1 = a.getWorldTransform()(pInA);
- supp->v2 = b.getWorldTransform()(qInB);
- supp->v = supp->v1 - supp->v2;
-}
-
-template <typename btConvexTemplate>
-static int btDiscoverPortal(const btConvexTemplate &a, const btConvexTemplate &b,
- const btMprCollisionDescription &colDesc,
- btMprSimplex_t *portal)
-{
- btVector3 dir, va, vb;
- float dot;
- int cont;
-
- // vertex 0 is center of portal
- btFindOrigin(a, b, colDesc, btMprSimplexPointW(portal, 0));
-
- // vertex 0 is center of portal
- btMprSimplexSetSize(portal, 1);
-
- btVector3 zero = btVector3(0, 0, 0);
- btVector3 *org = &zero;
-
- if (btMprVec3Eq(&btMprSimplexPoint(portal, 0)->v, org))
- {
- // Portal's center lies on origin (0,0,0) => we know that objects
- // intersect but we would need to know penetration info.
- // So move center little bit...
- btMprVec3Set(&va, FLT_EPSILON * 10.f, 0.f, 0.f);
- btMprVec3Add(&btMprSimplexPointW(portal, 0)->v, &va);
- }
-
- // vertex 1 = support in direction of origin
- btMprVec3Copy(&dir, &btMprSimplexPoint(portal, 0)->v);
- btMprVec3Scale(&dir, -1.f);
- btMprVec3Normalize(&dir);
-
- btMprSupport(a, b, colDesc, dir, btMprSimplexPointW(portal, 1));
-
- btMprSimplexSetSize(portal, 2);
-
- // test if origin isn't outside of v1
- dot = btMprVec3Dot(&btMprSimplexPoint(portal, 1)->v, &dir);
-
- if (btMprIsZero(dot) || dot < 0.f)
- return -1;
-
- // vertex 2
- btMprVec3Cross(&dir, &btMprSimplexPoint(portal, 0)->v,
- &btMprSimplexPoint(portal, 1)->v);
- if (btMprIsZero(btMprVec3Len2(&dir)))
- {
- if (btMprVec3Eq(&btMprSimplexPoint(portal, 1)->v, org))
- {
- // origin lies on v1
- return 1;
- }
- else
- {
- // origin lies on v0-v1 segment
- return 2;
- }
- }
-
- btMprVec3Normalize(&dir);
- btMprSupport(a, b, colDesc, dir, btMprSimplexPointW(portal, 2));
-
- dot = btMprVec3Dot(&btMprSimplexPoint(portal, 2)->v, &dir);
- if (btMprIsZero(dot) || dot < 0.f)
- return -1;
-
- btMprSimplexSetSize(portal, 3);
-
- // vertex 3 direction
- btMprVec3Sub2(&va, &btMprSimplexPoint(portal, 1)->v,
- &btMprSimplexPoint(portal, 0)->v);
- btMprVec3Sub2(&vb, &btMprSimplexPoint(portal, 2)->v,
- &btMprSimplexPoint(portal, 0)->v);
- btMprVec3Cross(&dir, &va, &vb);
- btMprVec3Normalize(&dir);
-
- // it is better to form portal faces to be oriented "outside" origin
- dot = btMprVec3Dot(&dir, &btMprSimplexPoint(portal, 0)->v);
- if (dot > 0.f)
- {
- btMprSimplexSwap(portal, 1, 2);
- btMprVec3Scale(&dir, -1.f);
- }
-
- while (btMprSimplexSize(portal) < 4)
- {
- btMprSupport(a, b, colDesc, dir, btMprSimplexPointW(portal, 3));
-
- dot = btMprVec3Dot(&btMprSimplexPoint(portal, 3)->v, &dir);
- if (btMprIsZero(dot) || dot < 0.f)
- return -1;
-
- cont = 0;
-
- // test if origin is outside (v1, v0, v3) - set v2 as v3 and
- // continue
- btMprVec3Cross(&va, &btMprSimplexPoint(portal, 1)->v,
- &btMprSimplexPoint(portal, 3)->v);
- dot = btMprVec3Dot(&va, &btMprSimplexPoint(portal, 0)->v);
- if (dot < 0.f && !btMprIsZero(dot))
- {
- btMprSimplexSet(portal, 2, btMprSimplexPoint(portal, 3));
- cont = 1;
- }
-
- if (!cont)
- {
- // test if origin is outside (v3, v0, v2) - set v1 as v3 and
- // continue
- btMprVec3Cross(&va, &btMprSimplexPoint(portal, 3)->v,
- &btMprSimplexPoint(portal, 2)->v);
- dot = btMprVec3Dot(&va, &btMprSimplexPoint(portal, 0)->v);
- if (dot < 0.f && !btMprIsZero(dot))
- {
- btMprSimplexSet(portal, 1, btMprSimplexPoint(portal, 3));
- cont = 1;
- }
- }
-
- if (cont)
- {
- btMprVec3Sub2(&va, &btMprSimplexPoint(portal, 1)->v,
- &btMprSimplexPoint(portal, 0)->v);
- btMprVec3Sub2(&vb, &btMprSimplexPoint(portal, 2)->v,
- &btMprSimplexPoint(portal, 0)->v);
- btMprVec3Cross(&dir, &va, &vb);
- btMprVec3Normalize(&dir);
- }
- else
- {
- btMprSimplexSetSize(portal, 4);
- }
- }
-
- return 0;
-}
-
-template <typename btConvexTemplate>
-static int btRefinePortal(const btConvexTemplate &a, const btConvexTemplate &b, const btMprCollisionDescription &colDesc,
- btMprSimplex_t *portal)
-{
- btVector3 dir;
- btMprSupport_t v4;
-
- for (int i = 0; i < BT_MPR_MAX_ITERATIONS; i++)
- //while (1)
- {
- // compute direction outside the portal (from v0 throught v1,v2,v3
- // face)
- btPortalDir(portal, &dir);
-
- // test if origin is inside the portal
- if (portalEncapsulesOrigin(portal, &dir))
- return 0;
-
- // get next support point
-
- btMprSupport(a, b, colDesc, dir, &v4);
-
- // test if v4 can expand portal to contain origin and if portal
- // expanding doesn't reach given tolerance
- if (!portalCanEncapsuleOrigin(portal, &v4, &dir) || portalReachTolerance(portal, &v4, &dir))
- {
- return -1;
- }
-
- // v1-v2-v3 triangle must be rearranged to face outside Minkowski
- // difference (direction from v0).
- btExpandPortal(portal, &v4);
- }
-
- return -1;
-}
-
-static void btFindPos(const btMprSimplex_t *portal, btVector3 *pos)
-{
- btVector3 zero = btVector3(0, 0, 0);
- btVector3 *origin = &zero;
-
- btVector3 dir;
- size_t i;
- float b[4], sum, inv;
- btVector3 vec, p1, p2;
-
- btPortalDir(portal, &dir);
-
- // use barycentric coordinates of tetrahedron to find origin
- btMprVec3Cross(&vec, &btMprSimplexPoint(portal, 1)->v,
- &btMprSimplexPoint(portal, 2)->v);
- b[0] = btMprVec3Dot(&vec, &btMprSimplexPoint(portal, 3)->v);
-
- btMprVec3Cross(&vec, &btMprSimplexPoint(portal, 3)->v,
- &btMprSimplexPoint(portal, 2)->v);
- b[1] = btMprVec3Dot(&vec, &btMprSimplexPoint(portal, 0)->v);
-
- btMprVec3Cross(&vec, &btMprSimplexPoint(portal, 0)->v,
- &btMprSimplexPoint(portal, 1)->v);
- b[2] = btMprVec3Dot(&vec, &btMprSimplexPoint(portal, 3)->v);
-
- btMprVec3Cross(&vec, &btMprSimplexPoint(portal, 2)->v,
- &btMprSimplexPoint(portal, 1)->v);
- b[3] = btMprVec3Dot(&vec, &btMprSimplexPoint(portal, 0)->v);
-
- sum = b[0] + b[1] + b[2] + b[3];
-
- if (btMprIsZero(sum) || sum < 0.f)
- {
- b[0] = 0.f;
-
- btMprVec3Cross(&vec, &btMprSimplexPoint(portal, 2)->v,
- &btMprSimplexPoint(portal, 3)->v);
- b[1] = btMprVec3Dot(&vec, &dir);
- btMprVec3Cross(&vec, &btMprSimplexPoint(portal, 3)->v,
- &btMprSimplexPoint(portal, 1)->v);
- b[2] = btMprVec3Dot(&vec, &dir);
- btMprVec3Cross(&vec, &btMprSimplexPoint(portal, 1)->v,
- &btMprSimplexPoint(portal, 2)->v);
- b[3] = btMprVec3Dot(&vec, &dir);
-
- sum = b[1] + b[2] + b[3];
- }
-
- inv = 1.f / sum;
-
- btMprVec3Copy(&p1, origin);
- btMprVec3Copy(&p2, origin);
- for (i = 0; i < 4; i++)
- {
- btMprVec3Copy(&vec, &btMprSimplexPoint(portal, i)->v1);
- btMprVec3Scale(&vec, b[i]);
- btMprVec3Add(&p1, &vec);
-
- btMprVec3Copy(&vec, &btMprSimplexPoint(portal, i)->v2);
- btMprVec3Scale(&vec, b[i]);
- btMprVec3Add(&p2, &vec);
- }
- btMprVec3Scale(&p1, inv);
- btMprVec3Scale(&p2, inv);
-#ifdef MPR_AVERAGE_CONTACT_POSITIONS
- btMprVec3Copy(pos, &p1);
- btMprVec3Add(pos, &p2);
- btMprVec3Scale(pos, 0.5);
-#else
- btMprVec3Copy(pos, &p2);
-#endif //MPR_AVERAGE_CONTACT_POSITIONS
-}
-
-inline float btMprVec3Dist2(const btVector3 *a, const btVector3 *b)
-{
- btVector3 ab;
- btMprVec3Sub2(&ab, a, b);
- return btMprVec3Len2(&ab);
-}
-
-inline float _btMprVec3PointSegmentDist2(const btVector3 *P,
- const btVector3 *x0,
- const btVector3 *b,
- btVector3 *witness)
-{
- // The computation comes from solving equation of segment:
- // S(t) = x0 + t.d
- // where - x0 is initial point of segment
- // - d is direction of segment from x0 (|d| > 0)
- // - t belongs to <0, 1> interval
- //
- // Than, distance from a segment to some point P can be expressed:
- // D(t) = |x0 + t.d - P|^2
- // which is distance from any point on segment. Minimization
- // of this function brings distance from P to segment.
- // Minimization of D(t) leads to simple quadratic equation that's
- // solving is straightforward.
- //
- // Bonus of this method is witness point for free.
-
- float dist, t;
- btVector3 d, a;
-
- // direction of segment
- btMprVec3Sub2(&d, b, x0);
-
- // precompute vector from P to x0
- btMprVec3Sub2(&a, x0, P);
-
- t = -1.f * btMprVec3Dot(&a, &d);
- t /= btMprVec3Len2(&d);
-
- if (t < 0.f || btMprIsZero(t))
- {
- dist = btMprVec3Dist2(x0, P);
- if (witness)
- btMprVec3Copy(witness, x0);
- }
- else if (t > 1.f || btMprEq(t, 1.f))
- {
- dist = btMprVec3Dist2(b, P);
- if (witness)
- btMprVec3Copy(witness, b);
- }
- else
- {
- if (witness)
- {
- btMprVec3Copy(witness, &d);
- btMprVec3Scale(witness, t);
- btMprVec3Add(witness, x0);
- dist = btMprVec3Dist2(witness, P);
- }
- else
- {
- // recycling variables
- btMprVec3Scale(&d, t);
- btMprVec3Add(&d, &a);
- dist = btMprVec3Len2(&d);
- }
- }
-
- return dist;
-}
-
-inline float btMprVec3PointTriDist2(const btVector3 *P,
- const btVector3 *x0, const btVector3 *B,
- const btVector3 *C,
- btVector3 *witness)
-{
- // Computation comes from analytic expression for triangle (x0, B, C)
- // T(s, t) = x0 + s.d1 + t.d2, where d1 = B - x0 and d2 = C - x0 and
- // Then equation for distance is:
- // D(s, t) = | T(s, t) - P |^2
- // This leads to minimization of quadratic function of two variables.
- // The solution from is taken only if s is between 0 and 1, t is
- // between 0 and 1 and t + s < 1, otherwise distance from segment is
- // computed.
-
- btVector3 d1, d2, a;
- float u, v, w, p, q, r;
- float s, t, dist, dist2;
- btVector3 witness2;
-
- btMprVec3Sub2(&d1, B, x0);
- btMprVec3Sub2(&d2, C, x0);
- btMprVec3Sub2(&a, x0, P);
-
- u = btMprVec3Dot(&a, &a);
- v = btMprVec3Dot(&d1, &d1);
- w = btMprVec3Dot(&d2, &d2);
- p = btMprVec3Dot(&a, &d1);
- q = btMprVec3Dot(&a, &d2);
- r = btMprVec3Dot(&d1, &d2);
-
- btScalar div = (w * v - r * r);
- if (btMprIsZero(div))
- {
- s = -1;
- }
- else
- {
- s = (q * r - w * p) / div;
- t = (-s * r - q) / w;
- }
-
- if ((btMprIsZero(s) || s > 0.f) && (btMprEq(s, 1.f) || s < 1.f) && (btMprIsZero(t) || t > 0.f) && (btMprEq(t, 1.f) || t < 1.f) && (btMprEq(t + s, 1.f) || t + s < 1.f))
- {
- if (witness)
- {
- btMprVec3Scale(&d1, s);
- btMprVec3Scale(&d2, t);
- btMprVec3Copy(witness, x0);
- btMprVec3Add(witness, &d1);
- btMprVec3Add(witness, &d2);
-
- dist = btMprVec3Dist2(witness, P);
- }
- else
- {
- dist = s * s * v;
- dist += t * t * w;
- dist += 2.f * s * t * r;
- dist += 2.f * s * p;
- dist += 2.f * t * q;
- dist += u;
- }
- }
- else
- {
- dist = _btMprVec3PointSegmentDist2(P, x0, B, witness);
-
- dist2 = _btMprVec3PointSegmentDist2(P, x0, C, &witness2);
- if (dist2 < dist)
- {
- dist = dist2;
- if (witness)
- btMprVec3Copy(witness, &witness2);
- }
-
- dist2 = _btMprVec3PointSegmentDist2(P, B, C, &witness2);
- if (dist2 < dist)
- {
- dist = dist2;
- if (witness)
- btMprVec3Copy(witness, &witness2);
- }
- }
-
- return dist;
-}
-
-template <typename btConvexTemplate>
-static void btFindPenetr(const btConvexTemplate &a, const btConvexTemplate &b,
- const btMprCollisionDescription &colDesc,
- btMprSimplex_t *portal,
- float *depth, btVector3 *pdir, btVector3 *pos)
-{
- btVector3 dir;
- btMprSupport_t v4;
- unsigned long iterations;
-
- btVector3 zero = btVector3(0, 0, 0);
- btVector3 *origin = &zero;
-
- iterations = 1UL;
- for (int i = 0; i < BT_MPR_MAX_ITERATIONS; i++)
- //while (1)
- {
- // compute portal direction and obtain next support point
- btPortalDir(portal, &dir);
-
- btMprSupport(a, b, colDesc, dir, &v4);
-
- // reached tolerance -> find penetration info
- if (portalReachTolerance(portal, &v4, &dir) || iterations == BT_MPR_MAX_ITERATIONS)
- {
- *depth = btMprVec3PointTriDist2(origin, &btMprSimplexPoint(portal, 1)->v, &btMprSimplexPoint(portal, 2)->v, &btMprSimplexPoint(portal, 3)->v, pdir);
- *depth = BT_MPR_SQRT(*depth);
-
- if (btMprIsZero((*pdir).x()) && btMprIsZero((*pdir).y()) && btMprIsZero((*pdir).z()))
- {
- *pdir = dir;
- }
- btMprVec3Normalize(pdir);
-
- // barycentric coordinates:
- btFindPos(portal, pos);
-
- return;
- }
-
- btExpandPortal(portal, &v4);
-
- iterations++;
- }
-}
-
-static void btFindPenetrTouch(btMprSimplex_t *portal, float *depth, btVector3 *dir, btVector3 *pos)
-{
- // Touching contact on portal's v1 - so depth is zero and direction
- // is unimportant and pos can be guessed
- *depth = 0.f;
- btVector3 zero = btVector3(0, 0, 0);
- btVector3 *origin = &zero;
-
- btMprVec3Copy(dir, origin);
-#ifdef MPR_AVERAGE_CONTACT_POSITIONS
- btMprVec3Copy(pos, &btMprSimplexPoint(portal, 1)->v1);
- btMprVec3Add(pos, &btMprSimplexPoint(portal, 1)->v2);
- btMprVec3Scale(pos, 0.5);
-#else
- btMprVec3Copy(pos, &btMprSimplexPoint(portal, 1)->v2);
-#endif
-}
-
-static void btFindPenetrSegment(btMprSimplex_t *portal,
- float *depth, btVector3 *dir, btVector3 *pos)
-{
- // Origin lies on v0-v1 segment.
- // Depth is distance to v1, direction also and position must be
- // computed
-#ifdef MPR_AVERAGE_CONTACT_POSITIONS
- btMprVec3Copy(pos, &btMprSimplexPoint(portal, 1)->v1);
- btMprVec3Add(pos, &btMprSimplexPoint(portal, 1)->v2);
- btMprVec3Scale(pos, 0.5f);
-#else
- btMprVec3Copy(pos, &btMprSimplexPoint(portal, 1)->v2);
-#endif //MPR_AVERAGE_CONTACT_POSITIONS
-
- btMprVec3Copy(dir, &btMprSimplexPoint(portal, 1)->v);
- *depth = BT_MPR_SQRT(btMprVec3Len2(dir));
- btMprVec3Normalize(dir);
-}
-
-template <typename btConvexTemplate>
-inline int btMprPenetration(const btConvexTemplate &a, const btConvexTemplate &b,
- const btMprCollisionDescription &colDesc,
- float *depthOut, btVector3 *dirOut, btVector3 *posOut)
-{
- btMprSimplex_t portal;
-
- // Phase 1: Portal discovery
- int result = btDiscoverPortal(a, b, colDesc, &portal);
-
- //sepAxis[pairIndex] = *pdir;//or -dir?
-
- switch (result)
- {
- case 0:
- {
- // Phase 2: Portal refinement
-
- result = btRefinePortal(a, b, colDesc, &portal);
- if (result < 0)
- return -1;
-
- // Phase 3. Penetration info
- btFindPenetr(a, b, colDesc, &portal, depthOut, dirOut, posOut);
-
- break;
- }
- case 1:
- {
- // Touching contact on portal's v1.
- btFindPenetrTouch(&portal, depthOut, dirOut, posOut);
- result = 0;
- break;
- }
- case 2:
- {
- btFindPenetrSegment(&portal, depthOut, dirOut, posOut);
- result = 0;
- break;
- }
- default:
- {
- //if (res < 0)
- //{
- // Origin isn't inside portal - no collision.
- result = -1;
- //}
- }
- };
-
- return result;
-};
-
-template <typename btConvexTemplate, typename btMprDistanceTemplate>
-inline int btComputeMprPenetration(const btConvexTemplate &a, const btConvexTemplate &b, const btMprCollisionDescription &colDesc, btMprDistanceTemplate *distInfo)
-{
- btVector3 dir, pos;
- float depth;
-
- int res = btMprPenetration(a, b, colDesc, &depth, &dir, &pos);
- if (res == 0)
- {
- distInfo->m_distance = -depth;
- distInfo->m_pointOnB = pos;
- distInfo->m_normalBtoA = -dir;
- distInfo->m_pointOnA = pos - distInfo->m_distance * dir;
- return 0;
- }
-
- return -1;
-}
-
-#endif //BT_MPR_PENETRATION_H
+ +/*** + * --------------------------------- + * Copyright (c)2012 Daniel Fiser <danfis@danfis.cz> + * + * This file was ported from mpr.c file, part of libccd. + * The Minkoski Portal Refinement implementation was ported + * to OpenCL by Erwin Coumans for the Bullet 3 Physics library. + * The original MPR idea and implementation is by Gary Snethen + * in XenoCollide, see http://github.com/erwincoumans/xenocollide + * + * Distributed under the OSI-approved BSD License (the "License"); + * see <http://www.opensource.org/licenses/bsd-license.php>. + * This software is distributed WITHOUT ANY WARRANTY; without even the + * implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + * See the License for more information. + */ + +///2014 Oct, Erwin Coumans, Use templates to avoid void* casts + +#ifndef BT_MPR_PENETRATION_H +#define BT_MPR_PENETRATION_H + +#define BT_DEBUG_MPR1 + +#include "LinearMath/btTransform.h" +#include "LinearMath/btAlignedObjectArray.h" + +//#define MPR_AVERAGE_CONTACT_POSITIONS + +struct btMprCollisionDescription +{ + btVector3 m_firstDir; + int m_maxGjkIterations; + btScalar m_maximumDistanceSquared; + btScalar m_gjkRelError2; + + btMprCollisionDescription() + : m_firstDir(0, 1, 0), + m_maxGjkIterations(1000), + m_maximumDistanceSquared(1e30f), + m_gjkRelError2(1.0e-6) + { + } + virtual ~btMprCollisionDescription() + { + } +}; + +struct btMprDistanceInfo +{ + btVector3 m_pointOnA; + btVector3 m_pointOnB; + btVector3 m_normalBtoA; + btScalar m_distance; +}; + +#ifdef __cplusplus +#define BT_MPR_SQRT sqrtf +#else +#define BT_MPR_SQRT sqrt +#endif +#define BT_MPR_FMIN(x, y) ((x) < (y) ? (x) : (y)) +#define BT_MPR_FABS fabs + +#define BT_MPR_TOLERANCE 1E-6f +#define BT_MPR_MAX_ITERATIONS 1000 + +struct _btMprSupport_t +{ + btVector3 v; //!< Support point in minkowski sum + btVector3 v1; //!< Support point in obj1 + btVector3 v2; //!< Support point in obj2 +}; +typedef struct _btMprSupport_t btMprSupport_t; + +struct _btMprSimplex_t +{ + btMprSupport_t ps[4]; + int last; //!< index of last added point +}; +typedef struct _btMprSimplex_t btMprSimplex_t; + +inline btMprSupport_t *btMprSimplexPointW(btMprSimplex_t *s, int idx) +{ + return &s->ps[idx]; +} + +inline void btMprSimplexSetSize(btMprSimplex_t *s, int size) +{ + s->last = size - 1; +} + +#ifdef DEBUG_MPR +inline void btPrintPortalVertex(_btMprSimplex_t *portal, int index) +{ + printf("portal[%d].v = %f,%f,%f, v1=%f,%f,%f, v2=%f,%f,%f\n", index, portal->ps[index].v.x(), portal->ps[index].v.y(), portal->ps[index].v.z(), + portal->ps[index].v1.x(), portal->ps[index].v1.y(), portal->ps[index].v1.z(), + portal->ps[index].v2.x(), portal->ps[index].v2.y(), portal->ps[index].v2.z()); +} +#endif //DEBUG_MPR + +inline int btMprSimplexSize(const btMprSimplex_t *s) +{ + return s->last + 1; +} + +inline const btMprSupport_t *btMprSimplexPoint(const btMprSimplex_t *s, int idx) +{ + // here is no check on boundaries + return &s->ps[idx]; +} + +inline void btMprSupportCopy(btMprSupport_t *d, const btMprSupport_t *s) +{ + *d = *s; +} + +inline void btMprSimplexSet(btMprSimplex_t *s, size_t pos, const btMprSupport_t *a) +{ + btMprSupportCopy(s->ps + pos, a); +} + +inline void btMprSimplexSwap(btMprSimplex_t *s, size_t pos1, size_t pos2) +{ + btMprSupport_t supp; + + btMprSupportCopy(&supp, &s->ps[pos1]); + btMprSupportCopy(&s->ps[pos1], &s->ps[pos2]); + btMprSupportCopy(&s->ps[pos2], &supp); +} + +inline int btMprIsZero(float val) +{ + return BT_MPR_FABS(val) < FLT_EPSILON; +} + +inline int btMprEq(float _a, float _b) +{ + float ab; + float a, b; + + ab = BT_MPR_FABS(_a - _b); + if (BT_MPR_FABS(ab) < FLT_EPSILON) + return 1; + + a = BT_MPR_FABS(_a); + b = BT_MPR_FABS(_b); + if (b > a) + { + return ab < FLT_EPSILON * b; + } + else + { + return ab < FLT_EPSILON * a; + } +} + +inline int btMprVec3Eq(const btVector3 *a, const btVector3 *b) +{ + return btMprEq((*a).x(), (*b).x()) && btMprEq((*a).y(), (*b).y()) && btMprEq((*a).z(), (*b).z()); +} + +template <typename btConvexTemplate> +inline void btFindOrigin(const btConvexTemplate &a, const btConvexTemplate &b, const btMprCollisionDescription &colDesc, btMprSupport_t *center) +{ + center->v1 = a.getObjectCenterInWorld(); + center->v2 = b.getObjectCenterInWorld(); + center->v = center->v1 - center->v2; +} + +inline void btMprVec3Set(btVector3 *v, float x, float y, float z) +{ + v->setValue(x, y, z); +} + +inline void btMprVec3Add(btVector3 *v, const btVector3 *w) +{ + *v += *w; +} + +inline void btMprVec3Copy(btVector3 *v, const btVector3 *w) +{ + *v = *w; +} + +inline void btMprVec3Scale(btVector3 *d, float k) +{ + *d *= k; +} + +inline float btMprVec3Dot(const btVector3 *a, const btVector3 *b) +{ + float dot; + + dot = btDot(*a, *b); + return dot; +} + +inline float btMprVec3Len2(const btVector3 *v) +{ + return btMprVec3Dot(v, v); +} + +inline void btMprVec3Normalize(btVector3 *d) +{ + float k = 1.f / BT_MPR_SQRT(btMprVec3Len2(d)); + btMprVec3Scale(d, k); +} + +inline void btMprVec3Cross(btVector3 *d, const btVector3 *a, const btVector3 *b) +{ + *d = btCross(*a, *b); +} + +inline void btMprVec3Sub2(btVector3 *d, const btVector3 *v, const btVector3 *w) +{ + *d = *v - *w; +} + +inline void btPortalDir(const btMprSimplex_t *portal, btVector3 *dir) +{ + btVector3 v2v1, v3v1; + + btMprVec3Sub2(&v2v1, &btMprSimplexPoint(portal, 2)->v, + &btMprSimplexPoint(portal, 1)->v); + btMprVec3Sub2(&v3v1, &btMprSimplexPoint(portal, 3)->v, + &btMprSimplexPoint(portal, 1)->v); + btMprVec3Cross(dir, &v2v1, &v3v1); + btMprVec3Normalize(dir); +} + +inline int portalEncapsulesOrigin(const btMprSimplex_t *portal, + const btVector3 *dir) +{ + float dot; + dot = btMprVec3Dot(dir, &btMprSimplexPoint(portal, 1)->v); + return btMprIsZero(dot) || dot > 0.f; +} + +inline int portalReachTolerance(const btMprSimplex_t *portal, + const btMprSupport_t *v4, + const btVector3 *dir) +{ + float dv1, dv2, dv3, dv4; + float dot1, dot2, dot3; + + // find the smallest dot product of dir and {v1-v4, v2-v4, v3-v4} + + dv1 = btMprVec3Dot(&btMprSimplexPoint(portal, 1)->v, dir); + dv2 = btMprVec3Dot(&btMprSimplexPoint(portal, 2)->v, dir); + dv3 = btMprVec3Dot(&btMprSimplexPoint(portal, 3)->v, dir); + dv4 = btMprVec3Dot(&v4->v, dir); + + dot1 = dv4 - dv1; + dot2 = dv4 - dv2; + dot3 = dv4 - dv3; + + dot1 = BT_MPR_FMIN(dot1, dot2); + dot1 = BT_MPR_FMIN(dot1, dot3); + + return btMprEq(dot1, BT_MPR_TOLERANCE) || dot1 < BT_MPR_TOLERANCE; +} + +inline int portalCanEncapsuleOrigin(const btMprSimplex_t *portal, + const btMprSupport_t *v4, + const btVector3 *dir) +{ + float dot; + dot = btMprVec3Dot(&v4->v, dir); + return btMprIsZero(dot) || dot > 0.f; +} + +inline void btExpandPortal(btMprSimplex_t *portal, + const btMprSupport_t *v4) +{ + float dot; + btVector3 v4v0; + + btMprVec3Cross(&v4v0, &v4->v, &btMprSimplexPoint(portal, 0)->v); + dot = btMprVec3Dot(&btMprSimplexPoint(portal, 1)->v, &v4v0); + if (dot > 0.f) + { + dot = btMprVec3Dot(&btMprSimplexPoint(portal, 2)->v, &v4v0); + if (dot > 0.f) + { + btMprSimplexSet(portal, 1, v4); + } + else + { + btMprSimplexSet(portal, 3, v4); + } + } + else + { + dot = btMprVec3Dot(&btMprSimplexPoint(portal, 3)->v, &v4v0); + if (dot > 0.f) + { + btMprSimplexSet(portal, 2, v4); + } + else + { + btMprSimplexSet(portal, 1, v4); + } + } +} +template <typename btConvexTemplate> +inline void btMprSupport(const btConvexTemplate &a, const btConvexTemplate &b, + const btMprCollisionDescription &colDesc, + const btVector3 &dir, btMprSupport_t *supp) +{ + btVector3 seperatingAxisInA = dir * a.getWorldTransform().getBasis(); + btVector3 seperatingAxisInB = -dir * b.getWorldTransform().getBasis(); + + btVector3 pInA = a.getLocalSupportWithMargin(seperatingAxisInA); + btVector3 qInB = b.getLocalSupportWithMargin(seperatingAxisInB); + + supp->v1 = a.getWorldTransform()(pInA); + supp->v2 = b.getWorldTransform()(qInB); + supp->v = supp->v1 - supp->v2; +} + +template <typename btConvexTemplate> +static int btDiscoverPortal(const btConvexTemplate &a, const btConvexTemplate &b, + const btMprCollisionDescription &colDesc, + btMprSimplex_t *portal) +{ + btVector3 dir, va, vb; + float dot; + int cont; + + // vertex 0 is center of portal + btFindOrigin(a, b, colDesc, btMprSimplexPointW(portal, 0)); + + // vertex 0 is center of portal + btMprSimplexSetSize(portal, 1); + + btVector3 zero = btVector3(0, 0, 0); + btVector3 *org = &zero; + + if (btMprVec3Eq(&btMprSimplexPoint(portal, 0)->v, org)) + { + // Portal's center lies on origin (0,0,0) => we know that objects + // intersect but we would need to know penetration info. + // So move center little bit... + btMprVec3Set(&va, FLT_EPSILON * 10.f, 0.f, 0.f); + btMprVec3Add(&btMprSimplexPointW(portal, 0)->v, &va); + } + + // vertex 1 = support in direction of origin + btMprVec3Copy(&dir, &btMprSimplexPoint(portal, 0)->v); + btMprVec3Scale(&dir, -1.f); + btMprVec3Normalize(&dir); + + btMprSupport(a, b, colDesc, dir, btMprSimplexPointW(portal, 1)); + + btMprSimplexSetSize(portal, 2); + + // test if origin isn't outside of v1 + dot = btMprVec3Dot(&btMprSimplexPoint(portal, 1)->v, &dir); + + if (btMprIsZero(dot) || dot < 0.f) + return -1; + + // vertex 2 + btMprVec3Cross(&dir, &btMprSimplexPoint(portal, 0)->v, + &btMprSimplexPoint(portal, 1)->v); + if (btMprIsZero(btMprVec3Len2(&dir))) + { + if (btMprVec3Eq(&btMprSimplexPoint(portal, 1)->v, org)) + { + // origin lies on v1 + return 1; + } + else + { + // origin lies on v0-v1 segment + return 2; + } + } + + btMprVec3Normalize(&dir); + btMprSupport(a, b, colDesc, dir, btMprSimplexPointW(portal, 2)); + + dot = btMprVec3Dot(&btMprSimplexPoint(portal, 2)->v, &dir); + if (btMprIsZero(dot) || dot < 0.f) + return -1; + + btMprSimplexSetSize(portal, 3); + + // vertex 3 direction + btMprVec3Sub2(&va, &btMprSimplexPoint(portal, 1)->v, + &btMprSimplexPoint(portal, 0)->v); + btMprVec3Sub2(&vb, &btMprSimplexPoint(portal, 2)->v, + &btMprSimplexPoint(portal, 0)->v); + btMprVec3Cross(&dir, &va, &vb); + btMprVec3Normalize(&dir); + + // it is better to form portal faces to be oriented "outside" origin + dot = btMprVec3Dot(&dir, &btMprSimplexPoint(portal, 0)->v); + if (dot > 0.f) + { + btMprSimplexSwap(portal, 1, 2); + btMprVec3Scale(&dir, -1.f); + } + + while (btMprSimplexSize(portal) < 4) + { + btMprSupport(a, b, colDesc, dir, btMprSimplexPointW(portal, 3)); + + dot = btMprVec3Dot(&btMprSimplexPoint(portal, 3)->v, &dir); + if (btMprIsZero(dot) || dot < 0.f) + return -1; + + cont = 0; + + // test if origin is outside (v1, v0, v3) - set v2 as v3 and + // continue + btMprVec3Cross(&va, &btMprSimplexPoint(portal, 1)->v, + &btMprSimplexPoint(portal, 3)->v); + dot = btMprVec3Dot(&va, &btMprSimplexPoint(portal, 0)->v); + if (dot < 0.f && !btMprIsZero(dot)) + { + btMprSimplexSet(portal, 2, btMprSimplexPoint(portal, 3)); + cont = 1; + } + + if (!cont) + { + // test if origin is outside (v3, v0, v2) - set v1 as v3 and + // continue + btMprVec3Cross(&va, &btMprSimplexPoint(portal, 3)->v, + &btMprSimplexPoint(portal, 2)->v); + dot = btMprVec3Dot(&va, &btMprSimplexPoint(portal, 0)->v); + if (dot < 0.f && !btMprIsZero(dot)) + { + btMprSimplexSet(portal, 1, btMprSimplexPoint(portal, 3)); + cont = 1; + } + } + + if (cont) + { + btMprVec3Sub2(&va, &btMprSimplexPoint(portal, 1)->v, + &btMprSimplexPoint(portal, 0)->v); + btMprVec3Sub2(&vb, &btMprSimplexPoint(portal, 2)->v, + &btMprSimplexPoint(portal, 0)->v); + btMprVec3Cross(&dir, &va, &vb); + btMprVec3Normalize(&dir); + } + else + { + btMprSimplexSetSize(portal, 4); + } + } + + return 0; +} + +template <typename btConvexTemplate> +static int btRefinePortal(const btConvexTemplate &a, const btConvexTemplate &b, const btMprCollisionDescription &colDesc, + btMprSimplex_t *portal) +{ + btVector3 dir; + btMprSupport_t v4; + + for (int i = 0; i < BT_MPR_MAX_ITERATIONS; i++) + //while (1) + { + // compute direction outside the portal (from v0 throught v1,v2,v3 + // face) + btPortalDir(portal, &dir); + + // test if origin is inside the portal + if (portalEncapsulesOrigin(portal, &dir)) + return 0; + + // get next support point + + btMprSupport(a, b, colDesc, dir, &v4); + + // test if v4 can expand portal to contain origin and if portal + // expanding doesn't reach given tolerance + if (!portalCanEncapsuleOrigin(portal, &v4, &dir) || portalReachTolerance(portal, &v4, &dir)) + { + return -1; + } + + // v1-v2-v3 triangle must be rearranged to face outside Minkowski + // difference (direction from v0). + btExpandPortal(portal, &v4); + } + + return -1; +} + +static void btFindPos(const btMprSimplex_t *portal, btVector3 *pos) +{ + btVector3 zero = btVector3(0, 0, 0); + btVector3 *origin = &zero; + + btVector3 dir; + size_t i; + float b[4], sum, inv; + btVector3 vec, p1, p2; + + btPortalDir(portal, &dir); + + // use barycentric coordinates of tetrahedron to find origin + btMprVec3Cross(&vec, &btMprSimplexPoint(portal, 1)->v, + &btMprSimplexPoint(portal, 2)->v); + b[0] = btMprVec3Dot(&vec, &btMprSimplexPoint(portal, 3)->v); + + btMprVec3Cross(&vec, &btMprSimplexPoint(portal, 3)->v, + &btMprSimplexPoint(portal, 2)->v); + b[1] = btMprVec3Dot(&vec, &btMprSimplexPoint(portal, 0)->v); + + btMprVec3Cross(&vec, &btMprSimplexPoint(portal, 0)->v, + &btMprSimplexPoint(portal, 1)->v); + b[2] = btMprVec3Dot(&vec, &btMprSimplexPoint(portal, 3)->v); + + btMprVec3Cross(&vec, &btMprSimplexPoint(portal, 2)->v, + &btMprSimplexPoint(portal, 1)->v); + b[3] = btMprVec3Dot(&vec, &btMprSimplexPoint(portal, 0)->v); + + sum = b[0] + b[1] + b[2] + b[3]; + + if (btMprIsZero(sum) || sum < 0.f) + { + b[0] = 0.f; + + btMprVec3Cross(&vec, &btMprSimplexPoint(portal, 2)->v, + &btMprSimplexPoint(portal, 3)->v); + b[1] = btMprVec3Dot(&vec, &dir); + btMprVec3Cross(&vec, &btMprSimplexPoint(portal, 3)->v, + &btMprSimplexPoint(portal, 1)->v); + b[2] = btMprVec3Dot(&vec, &dir); + btMprVec3Cross(&vec, &btMprSimplexPoint(portal, 1)->v, + &btMprSimplexPoint(portal, 2)->v); + b[3] = btMprVec3Dot(&vec, &dir); + + sum = b[1] + b[2] + b[3]; + } + + inv = 1.f / sum; + + btMprVec3Copy(&p1, origin); + btMprVec3Copy(&p2, origin); + for (i = 0; i < 4; i++) + { + btMprVec3Copy(&vec, &btMprSimplexPoint(portal, i)->v1); + btMprVec3Scale(&vec, b[i]); + btMprVec3Add(&p1, &vec); + + btMprVec3Copy(&vec, &btMprSimplexPoint(portal, i)->v2); + btMprVec3Scale(&vec, b[i]); + btMprVec3Add(&p2, &vec); + } + btMprVec3Scale(&p1, inv); + btMprVec3Scale(&p2, inv); +#ifdef MPR_AVERAGE_CONTACT_POSITIONS + btMprVec3Copy(pos, &p1); + btMprVec3Add(pos, &p2); + btMprVec3Scale(pos, 0.5); +#else + btMprVec3Copy(pos, &p2); +#endif //MPR_AVERAGE_CONTACT_POSITIONS +} + +inline float btMprVec3Dist2(const btVector3 *a, const btVector3 *b) +{ + btVector3 ab; + btMprVec3Sub2(&ab, a, b); + return btMprVec3Len2(&ab); +} + +inline float _btMprVec3PointSegmentDist2(const btVector3 *P, + const btVector3 *x0, + const btVector3 *b, + btVector3 *witness) +{ + // The computation comes from solving equation of segment: + // S(t) = x0 + t.d + // where - x0 is initial point of segment + // - d is direction of segment from x0 (|d| > 0) + // - t belongs to <0, 1> interval + // + // Than, distance from a segment to some point P can be expressed: + // D(t) = |x0 + t.d - P|^2 + // which is distance from any point on segment. Minimization + // of this function brings distance from P to segment. + // Minimization of D(t) leads to simple quadratic equation that's + // solving is straightforward. + // + // Bonus of this method is witness point for free. + + float dist, t; + btVector3 d, a; + + // direction of segment + btMprVec3Sub2(&d, b, x0); + + // precompute vector from P to x0 + btMprVec3Sub2(&a, x0, P); + + t = -1.f * btMprVec3Dot(&a, &d); + t /= btMprVec3Len2(&d); + + if (t < 0.f || btMprIsZero(t)) + { + dist = btMprVec3Dist2(x0, P); + if (witness) + btMprVec3Copy(witness, x0); + } + else if (t > 1.f || btMprEq(t, 1.f)) + { + dist = btMprVec3Dist2(b, P); + if (witness) + btMprVec3Copy(witness, b); + } + else + { + if (witness) + { + btMprVec3Copy(witness, &d); + btMprVec3Scale(witness, t); + btMprVec3Add(witness, x0); + dist = btMprVec3Dist2(witness, P); + } + else + { + // recycling variables + btMprVec3Scale(&d, t); + btMprVec3Add(&d, &a); + dist = btMprVec3Len2(&d); + } + } + + return dist; +} + +inline float btMprVec3PointTriDist2(const btVector3 *P, + const btVector3 *x0, const btVector3 *B, + const btVector3 *C, + btVector3 *witness) +{ + // Computation comes from analytic expression for triangle (x0, B, C) + // T(s, t) = x0 + s.d1 + t.d2, where d1 = B - x0 and d2 = C - x0 and + // Then equation for distance is: + // D(s, t) = | T(s, t) - P |^2 + // This leads to minimization of quadratic function of two variables. + // The solution from is taken only if s is between 0 and 1, t is + // between 0 and 1 and t + s < 1, otherwise distance from segment is + // computed. + + btVector3 d1, d2, a; + float u, v, w, p, q, r; + float s, t, dist, dist2; + btVector3 witness2; + + btMprVec3Sub2(&d1, B, x0); + btMprVec3Sub2(&d2, C, x0); + btMprVec3Sub2(&a, x0, P); + + u = btMprVec3Dot(&a, &a); + v = btMprVec3Dot(&d1, &d1); + w = btMprVec3Dot(&d2, &d2); + p = btMprVec3Dot(&a, &d1); + q = btMprVec3Dot(&a, &d2); + r = btMprVec3Dot(&d1, &d2); + + btScalar div = (w * v - r * r); + if (btMprIsZero(div)) + { + s = -1; + } + else + { + s = (q * r - w * p) / div; + t = (-s * r - q) / w; + } + + if ((btMprIsZero(s) || s > 0.f) && (btMprEq(s, 1.f) || s < 1.f) && (btMprIsZero(t) || t > 0.f) && (btMprEq(t, 1.f) || t < 1.f) && (btMprEq(t + s, 1.f) || t + s < 1.f)) + { + if (witness) + { + btMprVec3Scale(&d1, s); + btMprVec3Scale(&d2, t); + btMprVec3Copy(witness, x0); + btMprVec3Add(witness, &d1); + btMprVec3Add(witness, &d2); + + dist = btMprVec3Dist2(witness, P); + } + else + { + dist = s * s * v; + dist += t * t * w; + dist += 2.f * s * t * r; + dist += 2.f * s * p; + dist += 2.f * t * q; + dist += u; + } + } + else + { + dist = _btMprVec3PointSegmentDist2(P, x0, B, witness); + + dist2 = _btMprVec3PointSegmentDist2(P, x0, C, &witness2); + if (dist2 < dist) + { + dist = dist2; + if (witness) + btMprVec3Copy(witness, &witness2); + } + + dist2 = _btMprVec3PointSegmentDist2(P, B, C, &witness2); + if (dist2 < dist) + { + dist = dist2; + if (witness) + btMprVec3Copy(witness, &witness2); + } + } + + return dist; +} + +template <typename btConvexTemplate> +static void btFindPenetr(const btConvexTemplate &a, const btConvexTemplate &b, + const btMprCollisionDescription &colDesc, + btMprSimplex_t *portal, + float *depth, btVector3 *pdir, btVector3 *pos) +{ + btVector3 dir; + btMprSupport_t v4; + unsigned long iterations; + + btVector3 zero = btVector3(0, 0, 0); + btVector3 *origin = &zero; + + iterations = 1UL; + for (int i = 0; i < BT_MPR_MAX_ITERATIONS; i++) + //while (1) + { + // compute portal direction and obtain next support point + btPortalDir(portal, &dir); + + btMprSupport(a, b, colDesc, dir, &v4); + + // reached tolerance -> find penetration info + if (portalReachTolerance(portal, &v4, &dir) || iterations == BT_MPR_MAX_ITERATIONS) + { + *depth = btMprVec3PointTriDist2(origin, &btMprSimplexPoint(portal, 1)->v, &btMprSimplexPoint(portal, 2)->v, &btMprSimplexPoint(portal, 3)->v, pdir); + *depth = BT_MPR_SQRT(*depth); + + if (btMprIsZero((*pdir).x()) && btMprIsZero((*pdir).y()) && btMprIsZero((*pdir).z())) + { + *pdir = dir; + } + btMprVec3Normalize(pdir); + + // barycentric coordinates: + btFindPos(portal, pos); + + return; + } + + btExpandPortal(portal, &v4); + + iterations++; + } +} + +static void btFindPenetrTouch(btMprSimplex_t *portal, float *depth, btVector3 *dir, btVector3 *pos) +{ + // Touching contact on portal's v1 - so depth is zero and direction + // is unimportant and pos can be guessed + *depth = 0.f; + btVector3 zero = btVector3(0, 0, 0); + btVector3 *origin = &zero; + + btMprVec3Copy(dir, origin); +#ifdef MPR_AVERAGE_CONTACT_POSITIONS + btMprVec3Copy(pos, &btMprSimplexPoint(portal, 1)->v1); + btMprVec3Add(pos, &btMprSimplexPoint(portal, 1)->v2); + btMprVec3Scale(pos, 0.5); +#else + btMprVec3Copy(pos, &btMprSimplexPoint(portal, 1)->v2); +#endif +} + +static void btFindPenetrSegment(btMprSimplex_t *portal, + float *depth, btVector3 *dir, btVector3 *pos) +{ + // Origin lies on v0-v1 segment. + // Depth is distance to v1, direction also and position must be + // computed +#ifdef MPR_AVERAGE_CONTACT_POSITIONS + btMprVec3Copy(pos, &btMprSimplexPoint(portal, 1)->v1); + btMprVec3Add(pos, &btMprSimplexPoint(portal, 1)->v2); + btMprVec3Scale(pos, 0.5f); +#else + btMprVec3Copy(pos, &btMprSimplexPoint(portal, 1)->v2); +#endif //MPR_AVERAGE_CONTACT_POSITIONS + + btMprVec3Copy(dir, &btMprSimplexPoint(portal, 1)->v); + *depth = BT_MPR_SQRT(btMprVec3Len2(dir)); + btMprVec3Normalize(dir); +} + +template <typename btConvexTemplate> +inline int btMprPenetration(const btConvexTemplate &a, const btConvexTemplate &b, + const btMprCollisionDescription &colDesc, + float *depthOut, btVector3 *dirOut, btVector3 *posOut) +{ + btMprSimplex_t portal; + + // Phase 1: Portal discovery + int result = btDiscoverPortal(a, b, colDesc, &portal); + + //sepAxis[pairIndex] = *pdir;//or -dir? + + switch (result) + { + case 0: + { + // Phase 2: Portal refinement + + result = btRefinePortal(a, b, colDesc, &portal); + if (result < 0) + return -1; + + // Phase 3. Penetration info + btFindPenetr(a, b, colDesc, &portal, depthOut, dirOut, posOut); + + break; + } + case 1: + { + // Touching contact on portal's v1. + btFindPenetrTouch(&portal, depthOut, dirOut, posOut); + result = 0; + break; + } + case 2: + { + btFindPenetrSegment(&portal, depthOut, dirOut, posOut); + result = 0; + break; + } + default: + { + //if (res < 0) + //{ + // Origin isn't inside portal - no collision. + result = -1; + //} + } + }; + + return result; +}; + +template <typename btConvexTemplate, typename btMprDistanceTemplate> +inline int btComputeMprPenetration(const btConvexTemplate &a, const btConvexTemplate &b, const btMprCollisionDescription &colDesc, btMprDistanceTemplate *distInfo) +{ + btVector3 dir, pos; + float depth; + + int res = btMprPenetration(a, b, colDesc, &depth, &dir, &pos); + if (res == 0) + { + distInfo->m_distance = -depth; + distInfo->m_pointOnB = pos; + distInfo->m_normalBtoA = -dir; + distInfo->m_pointOnA = pos - distInfo->m_distance * dir; + return 0; + } + + return -1; +} + +#endif //BT_MPR_PENETRATION_H |