polycode-engine/include/BulletCollision/NarrowPhaseCollision/btMprPenetration.h

/***
 * ---------------------------------
 * Copyright (c)2012 Daniel Fiser <[email protected]>
 *
 *  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