GJK solver implementation from Bullet, adapted to C# for the HumanoidRetargeter project. It runs the Gilbert-Johnson-Keerthi iterative algorithm to compute distance or intersection (inside) between two convex shapes via a Minkowski difference support function, manages simplex vertices and projection helpers.
/* * C# / XNA port of Bullet (c) 2011 Mark Neale <[email protected]> * * Bullet Continuous Collision Detection and Physics Library * Copyright (c) 2003-2008 Erwin Coumans http://www.bulletphysics.com/ * * This software is provided 'as-is', without any express or implied warranty. * In no event will the authors be held liable for any damages arising from * the use of this software. * * Permission is granted to anyone to use this software for any purpose, * including commercial applications, and to alter it and redistribute it * freely, subject to the following restrictions: * * 1. The origin of this software must not be misrepresented; you must not * claim that you wrote the original software. If you use this software * in a product, an acknowledgment in the product documentation would be * appreciated but is not required. * 2. Altered source versions must be plainly marked as such, and must not be * misrepresented as being the original software. * 3. This notice may not be removed or altered from any source distribution. */ // Adapted from BulletXNA 73b2e1b9: namespace isolation, no XNA rendering, one type per file. #nullable disable #pragma warning disable CS1591, CS0168, CS0219, CS0414, CS0660, CS0661, CS0659, CS0693, CS8981, CS0108, CS0114 using System; using System.Collections.Generic; using HumanoidRetargeter.Core.Formats.Mmd.Bullet.LinearMath; using System.Diagnostics; namespace HumanoidRetargeter.Core.Formats.Mmd.Bullet.BulletCollision; internal class GJK : IDisposable { public GJK() { //Initialise(); } public void Initialise() { m_ray = IndexedVector3.Zero; m_nfree = 0; m_status = GJKStatus.Failed; m_current = 0; m_distance = 0f; for (int i = 0; i < m_simplices.Length; ++i) { if (m_simplices[i] == null) { m_simplices[i] = new sSimplex(); } } for (int i = 0; i < m_store.Length; ++i) { if (m_store[i] == null) { m_store[i] = new sSV(); } } } IndexedVector3[] lastw = new IndexedVector3[4]; public GJKStatus Evaluate(GjkEpaSolver2MinkowskiDiff shapearg, ref IndexedVector3 guess) { uint iterations = 0; float sqdist = 0f; float alpha = 0f; uint clastw = 0; /* Initialize solver */ m_free[0] = m_store[0]; m_free[1] = m_store[1]; m_free[2] = m_store[2]; m_free[3] = m_store[3]; m_nfree = 4; m_current = 0; m_status = GJKStatus.Valid; m_shape = shapearg; m_distance = 0f; /* Initialize simplex */ m_simplices[0].rank = 0; m_ray = guess; float sqrl = m_ray.LengthSquared(); IndexedVector3 temp = sqrl > 0 ? -m_ray : new IndexedVector3(1, 0, 0); AppendVertice(m_simplices[0], ref temp); m_simplices[0].p[0] = 1; m_ray = m_simplices[0].c[0].w; sqdist = sqrl; lastw[0] = lastw[1] = lastw[2] = lastw[3] = m_ray; /* Loop */ do { uint next = 1 - m_current; sSimplex cs = m_simplices[m_current]; sSimplex ns = m_simplices[next]; /* Check zero */ float rl = m_ray.Length(); if (rl < GjkEpaSolver2.GJK_MIN_DISTANCE) { /* Touching or inside */ m_status = GJKStatus.Inside; break; } /* Append new vertice in -'v' direction */ IndexedVector3 temp2 = -m_ray; AppendVertice(cs, ref temp2); IndexedVector3 w = cs.c[cs.rank - 1].w; bool found = false; for (int i = 0; i < 4; ++i) { if ((w - lastw[i]).LengthSquared() < GjkEpaSolver2.GJK_DUPLICATED_EPS) { found = true; break; } } if (found) { /* Return old simplex */ RemoveVertice(m_simplices[m_current]); break; } else { /* Update lastw */ lastw[clastw = (clastw + 1) & 3] = w; } /* Check for termination */ float omega = IndexedVector3.Dot(ref m_ray, ref w) / rl; alpha = Math.Max(omega, alpha); if (((rl - alpha) - (GjkEpaSolver2.GJK_ACCURARY * rl)) <= 0) { /* Return old simplex */ RemoveVertice(m_simplices[m_current]); break; } /* Reduce simplex */ IndexedVector4 weights = new IndexedVector4(); uint mask = 0; switch (cs.rank) { case 2: { sqdist = GJK.ProjectOrigin(ref cs.c[0].w, ref cs.c[1].w, ref weights, ref mask); break; } case 3: { sqdist = GJK.ProjectOrigin(ref cs.c[0].w, ref cs.c[1].w, ref cs.c[2].w, ref weights, ref mask); break; } case 4: { sqdist = GJK.ProjectOrigin(ref cs.c[0].w, ref cs.c[1].w, ref cs.c[2].w, ref cs.c[3].w, ref weights, ref mask); break; } } if (sqdist >= 0) { /* Valid */ ns.rank = 0; m_ray = IndexedVector3.Zero; m_current = next; for (uint i = 0, ni = cs.rank; i < ni; ++i) { if ((mask & (1 << (int)i)) != 0) { ns.c[ns.rank] = cs.c[i]; float weight = weights[(int)i]; ns.p[ns.rank++] = weight; m_ray += cs.c[i].w * weight; } else { m_free[m_nfree++] = cs.c[i]; } } if (mask == 15) { m_status = GJKStatus.Inside; } } else { /* Return old simplex */ RemoveVertice(m_simplices[m_current]); break; } m_status = ((++iterations) < GjkEpaSolver2.GJK_MAX_ITERATIONS) ? m_status : GJKStatus.Failed; } while (m_status == GJKStatus.Valid); m_simplex = m_simplices[m_current]; switch (m_status) { case GJKStatus.Valid: { m_distance = m_ray.Length(); break; } case GJKStatus.Inside: { m_distance = 0; break; } } return (m_status); } public bool EncloseOrigin() { switch (m_simplex.rank) { case 1: { for (int i = 0; i < 3; ++i) { IndexedVector3 axis = IndexedVector3.Zero; axis[i] = 1f; AppendVertice(m_simplex, ref axis); if (EncloseOrigin()) { return (true); } RemoveVertice(m_simplex); IndexedVector3 temp = -axis; AppendVertice(m_simplex, ref temp); if (EncloseOrigin()) { return (true); } RemoveVertice(m_simplex); } break; } case 2: { IndexedVector3 d = m_simplex.c[1].w - m_simplex.c[0].w; for (int i = 0; i < 3; ++i) { IndexedVector3 axis = IndexedVector3.Zero; axis[i] = 1f; IndexedVector3 p = IndexedVector3.Cross(d, axis); if (p.LengthSquared() > 0) { AppendVertice(m_simplex, ref p); if (EncloseOrigin()) { return (true); } RemoveVertice(m_simplex); IndexedVector3 temp = -p; AppendVertice(m_simplex, ref temp); if (EncloseOrigin()) { return (true); } RemoveVertice(m_simplex); } } break; } case 3: { IndexedVector3 n = IndexedVector3.Cross(m_simplex.c[1].w - m_simplex.c[0].w, m_simplex.c[2].w - m_simplex.c[0].w); if (n.LengthSquared() > 0) { AppendVertice(m_simplex, ref n); if (EncloseOrigin()) { return (true); } RemoveVertice(m_simplex); IndexedVector3 temp = -n; AppendVertice(m_simplex, ref temp); if (EncloseOrigin()) { return (true); } RemoveVertice(m_simplex); } break; } case 4: { if (Math.Abs(GJK.Det(m_simplex.c[0].w - m_simplex.c[3].w, m_simplex.c[1].w - m_simplex.c[3].w, m_simplex.c[2].w - m_simplex.c[3].w)) > 0) return (true); break; } default: { break; } } return (false); } public void GetSupport(ref IndexedVector3 d, ref sSV sv) { sv.d = d / d.Length(); sv.w = m_shape.Support(ref sv.d); } public void RemoveVertice(sSimplex simplex) { m_free[m_nfree++] = simplex.c[--simplex.rank]; } public void AppendVertice(sSimplex simplex, ref IndexedVector3 v) { simplex.p[simplex.rank] = 0; simplex.c[simplex.rank] = m_free[--m_nfree]; GetSupport(ref v, ref simplex.c[simplex.rank++]); } public static float Det(IndexedVector3 a, IndexedVector3 b, IndexedVector3 c) { return Det(ref a, ref b, ref c); } public static float Det(ref IndexedVector3 a, ref IndexedVector3 b, ref IndexedVector3 c) { return (a.Y * b.Z * c.X + a.Z * b.X * c.Y - a.X * b.Z * c.Y - a.Y * b.X * c.Z + a.X * b.Y * c.Z - a.Z * b.Y * c.X); } public static float ProjectOrigin(ref IndexedVector3 a, ref IndexedVector3 b, ref IndexedVector4 w, ref uint m) { IndexedVector3 d = b - a; float l = d.LengthSquared(); if (l > GjkEpaSolver2.GJK_SIMPLEX2_EPS) { float t = (l > 0f ? (-IndexedVector3.Dot(ref a, ref d) / l) : 0f); if (t >= 1) { w.X = 0f; w.Y = 1f; m = 2; return b.LengthSquared(); } else if (t <= 0) { w.X = 1f; w.Y = 0f; m = 1; return a.LengthSquared(); } else { w.X = 1 - (w.Y = t); m = 3; return (a + d * t).LengthSquared(); } } return (-1); } static uint[] imd3 = { 1, 2, 0 }; static IndexedVector3[] vt = new IndexedVector3[3]; static IndexedVector3[] dl = new IndexedVector3[3]; static bool inhere2 = false; public static float ProjectOrigin(ref IndexedVector3 a, ref IndexedVector3 b, ref IndexedVector3 c, ref IndexedVector4 w, ref uint m) { { } inhere2 = true; vt[0] = a; vt[1] = b; vt[2] = c; dl[0] = a - b; dl[1] = b - c; dl[2] = c - a; IndexedVector3 n = IndexedVector3.Cross(dl[0], dl[1]); float l = n.LengthSquared(); if (l > GjkEpaSolver2.GJK_SIMPLEX3_EPS) { float mindist = -1f; IndexedVector4 subw = new IndexedVector4(); uint subm = 0; for (int i = 0; i < 3; ++i) { if (IndexedVector3.Dot(vt[i], IndexedVector3.Cross(dl[i], n)) > 0) { uint j = imd3[i]; float subd = GJK.ProjectOrigin(ref vt[i], ref vt[j], ref subw, ref subm); if ((mindist < 0) || (subd < mindist)) { mindist = subd; m = (uint)((((subm & 1) != 0) ? 1 << i : 0) + (((subm & 2) != 0) ? 1 << (int)j : 0)); w[(int)i] = subw.X; w[(int)j] = subw.Y; w[(int)imd3[j]] = 0f; } } } if (mindist < 0) { float d = IndexedVector3.Dot(ref a, ref n); float s = (float)Math.Sqrt(l); IndexedVector3 p = n * (d / l); mindist = p.LengthSquared(); m = 7; w.X = (IndexedVector3.Cross(dl[1], b - p)).Length() / s; w.Y = (IndexedVector3.Cross(dl[2], c - p)).Length() / s; w.Z = 1 - (w.X + w.Y); } inhere2 = false; return (mindist); } inhere2 = false; return (-1); } private static uint[] imd3a = { 1, 2, 0 }; private static IndexedVector3[] vta = new IndexedVector3[4]; private static IndexedVector3[] dla = new IndexedVector3[3]; private static bool inhere1 = false; public static float ProjectOrigin(ref IndexedVector3 a, ref IndexedVector3 b, ref IndexedVector3 c, ref IndexedVector3 d, ref IndexedVector4 w, ref uint m) { { } inhere1 = true; vta[0] = a; vta[1] = b; vta[2] = c; vta[3] = d; dla[0] = a - d; dla[1] = b - d; dla[2] = c - d; float vl = Det(dl[0], dl[1], dl[2]); bool ng = (vl * IndexedVector3.Dot(a, IndexedVector3.Cross(b - c, a - b))) <= 0; if (ng && (Math.Abs(vl) > GjkEpaSolver2.GJK_SIMPLEX4_EPS)) { float mindist = -1; IndexedVector4 subw = new IndexedVector4(); uint subm = 0; for (int i = 0; i < 3; ++i) { uint j = imd3[i]; float s = vl * IndexedVector3.Dot(d, IndexedVector3.Cross(dl[i], dl[j])); if (s > 0) { float subd = GJK.ProjectOrigin(ref vt[i], ref vt[j], ref d, ref subw, ref subm); if ((mindist < 0) || (subd < mindist)) { mindist = subd; m = (uint)((((subm & 1) != 0) ? 1 << i : 0) + (((subm & 2) != 0) ? 1 << (int)j : 0) + (((subm & 4) != 0) ? 8 : 0)); w[(int)i] = subw.X; w[(int)j] = subw.Y; w[(int)imd3[j]] = 0f; w.W = subw.Z; } } } if (mindist < 0) { mindist = 0; m = 15; w.X = Det(c, b, d) / vl; w.Y = Det(a, c, d) / vl; w.Z = Det(b, a, d) / vl; w.W = 1 - (w.X + w.Y + w.Z); } inhere1 = false; return (mindist); } inhere1 = false; return (-1); } public void Dispose() { BulletGlobals.GJKPool.Free(this); } public GjkEpaSolver2MinkowskiDiff m_shape; public IndexedVector3 m_ray; public float m_distance; public sSimplex[] m_simplices = new sSimplex[2]; public sSV[] m_store = new sSV[4]; public sSV[] m_free = new sSV[4]; public uint m_nfree; public uint m_current; public sSimplex m_simplex = new sSimplex(); public GJKStatus m_status; }