Você não pode selecionar mais de 25 tópicos Os tópicos devem começar com uma letra ou um número, podem incluir traços ('-') e podem ter até 35 caracteres.
 
 
 

921 linhas
26 KiB

  1. /***
  2. * ---------------------------------
  3. * Copyright (c)2012 Daniel Fiser <danfis@danfis.cz>
  4. *
  5. * This file was ported from mpr.c file, part of libccd.
  6. * The Minkoski Portal Refinement implementation was ported
  7. * to OpenCL by Erwin Coumans for the Bullet 3 Physics library.
  8. * at http://github.com/erwincoumans/bullet3
  9. *
  10. * Distributed under the OSI-approved BSD License (the "License");
  11. * see <http://www.opensource.org/licenses/bsd-license.php>.
  12. * This software is distributed WITHOUT ANY WARRANTY; without even the
  13. * implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
  14. * See the License for more information.
  15. */
  16. #ifndef B3_MPR_PENETRATION_H
  17. #define B3_MPR_PENETRATION_H
  18. #include "Bullet3Common/shared/b3PlatformDefinitions.h"
  19. #include "Bullet3Common/shared/b3Float4.h"
  20. #include "Bullet3Collision/NarrowPhaseCollision/shared/b3RigidBodyData.h"
  21. #include "Bullet3Collision/NarrowPhaseCollision/shared/b3ConvexPolyhedronData.h"
  22. #include "Bullet3Collision/NarrowPhaseCollision/shared/b3Collidable.h"
  23. #ifdef __cplusplus
  24. #define B3_MPR_SQRT sqrtf
  25. #else
  26. #define B3_MPR_SQRT sqrt
  27. #endif
  28. #define B3_MPR_FMIN(x, y) ((x) < (y) ? (x) : (y))
  29. #define B3_MPR_FABS fabs
  30. #define B3_MPR_TOLERANCE 1E-6f
  31. #define B3_MPR_MAX_ITERATIONS 1000
  32. struct _b3MprSupport_t
  33. {
  34. b3Float4 v; //!< Support point in minkowski sum
  35. b3Float4 v1; //!< Support point in obj1
  36. b3Float4 v2; //!< Support point in obj2
  37. };
  38. typedef struct _b3MprSupport_t b3MprSupport_t;
  39. struct _b3MprSimplex_t
  40. {
  41. b3MprSupport_t ps[4];
  42. int last; //!< index of last added point
  43. };
  44. typedef struct _b3MprSimplex_t b3MprSimplex_t;
  45. inline b3MprSupport_t* b3MprSimplexPointW(b3MprSimplex_t *s, int idx)
  46. {
  47. return &s->ps[idx];
  48. }
  49. inline void b3MprSimplexSetSize(b3MprSimplex_t *s, int size)
  50. {
  51. s->last = size - 1;
  52. }
  53. inline int b3MprSimplexSize(const b3MprSimplex_t *s)
  54. {
  55. return s->last + 1;
  56. }
  57. inline const b3MprSupport_t* b3MprSimplexPoint(const b3MprSimplex_t* s, int idx)
  58. {
  59. // here is no check on boundaries
  60. return &s->ps[idx];
  61. }
  62. inline void b3MprSupportCopy(b3MprSupport_t *d, const b3MprSupport_t *s)
  63. {
  64. *d = *s;
  65. }
  66. inline void b3MprSimplexSet(b3MprSimplex_t *s, size_t pos, const b3MprSupport_t *a)
  67. {
  68. b3MprSupportCopy(s->ps + pos, a);
  69. }
  70. inline void b3MprSimplexSwap(b3MprSimplex_t *s, size_t pos1, size_t pos2)
  71. {
  72. b3MprSupport_t supp;
  73. b3MprSupportCopy(&supp, &s->ps[pos1]);
  74. b3MprSupportCopy(&s->ps[pos1], &s->ps[pos2]);
  75. b3MprSupportCopy(&s->ps[pos2], &supp);
  76. }
  77. inline int b3MprIsZero(float val)
  78. {
  79. return B3_MPR_FABS(val) < FLT_EPSILON;
  80. }
  81. inline int b3MprEq(float _a, float _b)
  82. {
  83. float ab;
  84. float a, b;
  85. ab = B3_MPR_FABS(_a - _b);
  86. if (B3_MPR_FABS(ab) < FLT_EPSILON)
  87. return 1;
  88. a = B3_MPR_FABS(_a);
  89. b = B3_MPR_FABS(_b);
  90. if (b > a){
  91. return ab < FLT_EPSILON * b;
  92. }else{
  93. return ab < FLT_EPSILON * a;
  94. }
  95. }
  96. inline int b3MprVec3Eq(const b3Float4* a, const b3Float4 *b)
  97. {
  98. return b3MprEq((*a).x, (*b).x)
  99. && b3MprEq((*a).y, (*b).y)
  100. && b3MprEq((*a).z, (*b).z);
  101. }
  102. inline b3Float4 b3LocalGetSupportVertex(b3Float4ConstArg supportVec,__global const b3ConvexPolyhedronData_t* hull, b3ConstArray(b3Float4) verticesA)
  103. {
  104. b3Float4 supVec = b3MakeFloat4(0,0,0,0);
  105. float maxDot = -B3_LARGE_FLOAT;
  106. if( 0 < hull->m_numVertices )
  107. {
  108. const b3Float4 scaled = supportVec;
  109. int index = b3MaxDot(scaled, &verticesA[hull->m_vertexOffset], hull->m_numVertices, &maxDot);
  110. return verticesA[hull->m_vertexOffset+index];
  111. }
  112. return supVec;
  113. }
  114. B3_STATIC void b3MprConvexSupport(int pairIndex,int bodyIndex, b3ConstArray(b3RigidBodyData_t) cpuBodyBuf,
  115. b3ConstArray(b3ConvexPolyhedronData_t) cpuConvexData,
  116. b3ConstArray(b3Collidable_t) cpuCollidables,
  117. b3ConstArray(b3Float4) cpuVertices,
  118. __global b3Float4* sepAxis,
  119. const b3Float4* _dir, b3Float4* outp, int logme)
  120. {
  121. //dir is in worldspace, move to local space
  122. b3Float4 pos = cpuBodyBuf[bodyIndex].m_pos;
  123. b3Quat orn = cpuBodyBuf[bodyIndex].m_quat;
  124. b3Float4 dir = b3MakeFloat4((*_dir).x,(*_dir).y,(*_dir).z,0.f);
  125. const b3Float4 localDir = b3QuatRotate(b3QuatInverse(orn),dir);
  126. //find local support vertex
  127. int colIndex = cpuBodyBuf[bodyIndex].m_collidableIdx;
  128. b3Assert(cpuCollidables[colIndex].m_shapeType==SHAPE_CONVEX_HULL);
  129. __global const b3ConvexPolyhedronData_t* hull = &cpuConvexData[cpuCollidables[colIndex].m_shapeIndex];
  130. b3Float4 pInA;
  131. if (logme)
  132. {
  133. b3Float4 supVec = b3MakeFloat4(0,0,0,0);
  134. float maxDot = -B3_LARGE_FLOAT;
  135. if( 0 < hull->m_numVertices )
  136. {
  137. const b3Float4 scaled = localDir;
  138. int index = b3MaxDot(scaled, &cpuVertices[hull->m_vertexOffset], hull->m_numVertices, &maxDot);
  139. pInA = cpuVertices[hull->m_vertexOffset+index];
  140. }
  141. } else
  142. {
  143. pInA = b3LocalGetSupportVertex(localDir,hull,cpuVertices);
  144. }
  145. //move vertex to world space
  146. *outp = b3TransformPoint(pInA,pos,orn);
  147. }
  148. inline void b3MprSupport(int pairIndex,int bodyIndexA, int bodyIndexB, b3ConstArray(b3RigidBodyData_t) cpuBodyBuf,
  149. b3ConstArray(b3ConvexPolyhedronData_t) cpuConvexData,
  150. b3ConstArray(b3Collidable_t) cpuCollidables,
  151. b3ConstArray(b3Float4) cpuVertices,
  152. __global b3Float4* sepAxis,
  153. const b3Float4* _dir, b3MprSupport_t *supp)
  154. {
  155. b3Float4 dir;
  156. dir = *_dir;
  157. b3MprConvexSupport(pairIndex,bodyIndexA,cpuBodyBuf,cpuConvexData,cpuCollidables,cpuVertices,sepAxis,&dir, &supp->v1,0);
  158. dir = *_dir*-1.f;
  159. b3MprConvexSupport(pairIndex,bodyIndexB,cpuBodyBuf,cpuConvexData,cpuCollidables,cpuVertices,sepAxis,&dir, &supp->v2,0);
  160. supp->v = supp->v1 - supp->v2;
  161. }
  162. inline void b3FindOrigin(int bodyIndexA, int bodyIndexB, b3ConstArray(b3RigidBodyData_t) cpuBodyBuf, b3MprSupport_t *center)
  163. {
  164. center->v1 = cpuBodyBuf[bodyIndexA].m_pos;
  165. center->v2 = cpuBodyBuf[bodyIndexB].m_pos;
  166. center->v = center->v1 - center->v2;
  167. }
  168. inline void b3MprVec3Set(b3Float4 *v, float x, float y, float z)
  169. {
  170. (*v).x = x;
  171. (*v).y = y;
  172. (*v).z = z;
  173. (*v).w = 0.f;
  174. }
  175. inline void b3MprVec3Add(b3Float4 *v, const b3Float4 *w)
  176. {
  177. (*v).x += (*w).x;
  178. (*v).y += (*w).y;
  179. (*v).z += (*w).z;
  180. }
  181. inline void b3MprVec3Copy(b3Float4 *v, const b3Float4 *w)
  182. {
  183. *v = *w;
  184. }
  185. inline void b3MprVec3Scale(b3Float4 *d, float k)
  186. {
  187. *d *= k;
  188. }
  189. inline float b3MprVec3Dot(const b3Float4 *a, const b3Float4 *b)
  190. {
  191. float dot;
  192. dot = b3Dot3F4(*a,*b);
  193. return dot;
  194. }
  195. inline float b3MprVec3Len2(const b3Float4 *v)
  196. {
  197. return b3MprVec3Dot(v, v);
  198. }
  199. inline void b3MprVec3Normalize(b3Float4 *d)
  200. {
  201. float k = 1.f / B3_MPR_SQRT(b3MprVec3Len2(d));
  202. b3MprVec3Scale(d, k);
  203. }
  204. inline void b3MprVec3Cross(b3Float4 *d, const b3Float4 *a, const b3Float4 *b)
  205. {
  206. *d = b3Cross3(*a,*b);
  207. }
  208. inline void b3MprVec3Sub2(b3Float4 *d, const b3Float4 *v, const b3Float4 *w)
  209. {
  210. *d = *v - *w;
  211. }
  212. inline void b3PortalDir(const b3MprSimplex_t *portal, b3Float4 *dir)
  213. {
  214. b3Float4 v2v1, v3v1;
  215. b3MprVec3Sub2(&v2v1, &b3MprSimplexPoint(portal, 2)->v,
  216. &b3MprSimplexPoint(portal, 1)->v);
  217. b3MprVec3Sub2(&v3v1, &b3MprSimplexPoint(portal, 3)->v,
  218. &b3MprSimplexPoint(portal, 1)->v);
  219. b3MprVec3Cross(dir, &v2v1, &v3v1);
  220. b3MprVec3Normalize(dir);
  221. }
  222. inline int portalEncapsulesOrigin(const b3MprSimplex_t *portal,
  223. const b3Float4 *dir)
  224. {
  225. float dot;
  226. dot = b3MprVec3Dot(dir, &b3MprSimplexPoint(portal, 1)->v);
  227. return b3MprIsZero(dot) || dot > 0.f;
  228. }
  229. inline int portalReachTolerance(const b3MprSimplex_t *portal,
  230. const b3MprSupport_t *v4,
  231. const b3Float4 *dir)
  232. {
  233. float dv1, dv2, dv3, dv4;
  234. float dot1, dot2, dot3;
  235. // find the smallest dot product of dir and {v1-v4, v2-v4, v3-v4}
  236. dv1 = b3MprVec3Dot(&b3MprSimplexPoint(portal, 1)->v, dir);
  237. dv2 = b3MprVec3Dot(&b3MprSimplexPoint(portal, 2)->v, dir);
  238. dv3 = b3MprVec3Dot(&b3MprSimplexPoint(portal, 3)->v, dir);
  239. dv4 = b3MprVec3Dot(&v4->v, dir);
  240. dot1 = dv4 - dv1;
  241. dot2 = dv4 - dv2;
  242. dot3 = dv4 - dv3;
  243. dot1 = B3_MPR_FMIN(dot1, dot2);
  244. dot1 = B3_MPR_FMIN(dot1, dot3);
  245. return b3MprEq(dot1, B3_MPR_TOLERANCE) || dot1 < B3_MPR_TOLERANCE;
  246. }
  247. inline int portalCanEncapsuleOrigin(const b3MprSimplex_t *portal,
  248. const b3MprSupport_t *v4,
  249. const b3Float4 *dir)
  250. {
  251. float dot;
  252. dot = b3MprVec3Dot(&v4->v, dir);
  253. return b3MprIsZero(dot) || dot > 0.f;
  254. }
  255. inline void b3ExpandPortal(b3MprSimplex_t *portal,
  256. const b3MprSupport_t *v4)
  257. {
  258. float dot;
  259. b3Float4 v4v0;
  260. b3MprVec3Cross(&v4v0, &v4->v, &b3MprSimplexPoint(portal, 0)->v);
  261. dot = b3MprVec3Dot(&b3MprSimplexPoint(portal, 1)->v, &v4v0);
  262. if (dot > 0.f){
  263. dot = b3MprVec3Dot(&b3MprSimplexPoint(portal, 2)->v, &v4v0);
  264. if (dot > 0.f){
  265. b3MprSimplexSet(portal, 1, v4);
  266. }else{
  267. b3MprSimplexSet(portal, 3, v4);
  268. }
  269. }else{
  270. dot = b3MprVec3Dot(&b3MprSimplexPoint(portal, 3)->v, &v4v0);
  271. if (dot > 0.f){
  272. b3MprSimplexSet(portal, 2, v4);
  273. }else{
  274. b3MprSimplexSet(portal, 1, v4);
  275. }
  276. }
  277. }
  278. B3_STATIC int b3DiscoverPortal(int pairIndex, int bodyIndexA, int bodyIndexB, b3ConstArray(b3RigidBodyData_t) cpuBodyBuf,
  279. b3ConstArray(b3ConvexPolyhedronData_t) cpuConvexData,
  280. b3ConstArray(b3Collidable_t) cpuCollidables,
  281. b3ConstArray(b3Float4) cpuVertices,
  282. __global b3Float4* sepAxis,
  283. __global int* hasSepAxis,
  284. b3MprSimplex_t *portal)
  285. {
  286. b3Float4 dir, va, vb;
  287. float dot;
  288. int cont;
  289. // vertex 0 is center of portal
  290. b3FindOrigin(bodyIndexA,bodyIndexB,cpuBodyBuf, b3MprSimplexPointW(portal, 0));
  291. // vertex 0 is center of portal
  292. b3MprSimplexSetSize(portal, 1);
  293. b3Float4 zero = b3MakeFloat4(0,0,0,0);
  294. b3Float4* b3mpr_vec3_origin = &zero;
  295. if (b3MprVec3Eq(&b3MprSimplexPoint(portal, 0)->v, b3mpr_vec3_origin)){
  296. // Portal's center lies on origin (0,0,0) => we know that objects
  297. // intersect but we would need to know penetration info.
  298. // So move center little bit...
  299. b3MprVec3Set(&va, FLT_EPSILON * 10.f, 0.f, 0.f);
  300. b3MprVec3Add(&b3MprSimplexPointW(portal, 0)->v, &va);
  301. }
  302. // vertex 1 = support in direction of origin
  303. b3MprVec3Copy(&dir, &b3MprSimplexPoint(portal, 0)->v);
  304. b3MprVec3Scale(&dir, -1.f);
  305. b3MprVec3Normalize(&dir);
  306. b3MprSupport(pairIndex,bodyIndexA,bodyIndexB,cpuBodyBuf,cpuConvexData,cpuCollidables,cpuVertices, sepAxis,&dir, b3MprSimplexPointW(portal, 1));
  307. b3MprSimplexSetSize(portal, 2);
  308. // test if origin isn't outside of v1
  309. dot = b3MprVec3Dot(&b3MprSimplexPoint(portal, 1)->v, &dir);
  310. if (b3MprIsZero(dot) || dot < 0.f)
  311. return -1;
  312. // vertex 2
  313. b3MprVec3Cross(&dir, &b3MprSimplexPoint(portal, 0)->v,
  314. &b3MprSimplexPoint(portal, 1)->v);
  315. if (b3MprIsZero(b3MprVec3Len2(&dir))){
  316. if (b3MprVec3Eq(&b3MprSimplexPoint(portal, 1)->v, b3mpr_vec3_origin)){
  317. // origin lies on v1
  318. return 1;
  319. }else{
  320. // origin lies on v0-v1 segment
  321. return 2;
  322. }
  323. }
  324. b3MprVec3Normalize(&dir);
  325. b3MprSupport(pairIndex,bodyIndexA,bodyIndexB,cpuBodyBuf,cpuConvexData,cpuCollidables,cpuVertices, sepAxis,&dir, b3MprSimplexPointW(portal, 2));
  326. dot = b3MprVec3Dot(&b3MprSimplexPoint(portal, 2)->v, &dir);
  327. if (b3MprIsZero(dot) || dot < 0.f)
  328. return -1;
  329. b3MprSimplexSetSize(portal, 3);
  330. // vertex 3 direction
  331. b3MprVec3Sub2(&va, &b3MprSimplexPoint(portal, 1)->v,
  332. &b3MprSimplexPoint(portal, 0)->v);
  333. b3MprVec3Sub2(&vb, &b3MprSimplexPoint(portal, 2)->v,
  334. &b3MprSimplexPoint(portal, 0)->v);
  335. b3MprVec3Cross(&dir, &va, &vb);
  336. b3MprVec3Normalize(&dir);
  337. // it is better to form portal faces to be oriented "outside" origin
  338. dot = b3MprVec3Dot(&dir, &b3MprSimplexPoint(portal, 0)->v);
  339. if (dot > 0.f){
  340. b3MprSimplexSwap(portal, 1, 2);
  341. b3MprVec3Scale(&dir, -1.f);
  342. }
  343. while (b3MprSimplexSize(portal) < 4){
  344. b3MprSupport(pairIndex,bodyIndexA,bodyIndexB,cpuBodyBuf,cpuConvexData,cpuCollidables,cpuVertices, sepAxis,&dir, b3MprSimplexPointW(portal, 3));
  345. dot = b3MprVec3Dot(&b3MprSimplexPoint(portal, 3)->v, &dir);
  346. if (b3MprIsZero(dot) || dot < 0.f)
  347. return -1;
  348. cont = 0;
  349. // test if origin is outside (v1, v0, v3) - set v2 as v3 and
  350. // continue
  351. b3MprVec3Cross(&va, &b3MprSimplexPoint(portal, 1)->v,
  352. &b3MprSimplexPoint(portal, 3)->v);
  353. dot = b3MprVec3Dot(&va, &b3MprSimplexPoint(portal, 0)->v);
  354. if (dot < 0.f && !b3MprIsZero(dot)){
  355. b3MprSimplexSet(portal, 2, b3MprSimplexPoint(portal, 3));
  356. cont = 1;
  357. }
  358. if (!cont){
  359. // test if origin is outside (v3, v0, v2) - set v1 as v3 and
  360. // continue
  361. b3MprVec3Cross(&va, &b3MprSimplexPoint(portal, 3)->v,
  362. &b3MprSimplexPoint(portal, 2)->v);
  363. dot = b3MprVec3Dot(&va, &b3MprSimplexPoint(portal, 0)->v);
  364. if (dot < 0.f && !b3MprIsZero(dot)){
  365. b3MprSimplexSet(portal, 1, b3MprSimplexPoint(portal, 3));
  366. cont = 1;
  367. }
  368. }
  369. if (cont){
  370. b3MprVec3Sub2(&va, &b3MprSimplexPoint(portal, 1)->v,
  371. &b3MprSimplexPoint(portal, 0)->v);
  372. b3MprVec3Sub2(&vb, &b3MprSimplexPoint(portal, 2)->v,
  373. &b3MprSimplexPoint(portal, 0)->v);
  374. b3MprVec3Cross(&dir, &va, &vb);
  375. b3MprVec3Normalize(&dir);
  376. }else{
  377. b3MprSimplexSetSize(portal, 4);
  378. }
  379. }
  380. return 0;
  381. }
  382. B3_STATIC int b3RefinePortal(int pairIndex,int bodyIndexA, int bodyIndexB, b3ConstArray(b3RigidBodyData_t) cpuBodyBuf,
  383. b3ConstArray(b3ConvexPolyhedronData_t) cpuConvexData,
  384. b3ConstArray(b3Collidable_t) cpuCollidables,
  385. b3ConstArray(b3Float4) cpuVertices,
  386. __global b3Float4* sepAxis,
  387. b3MprSimplex_t *portal)
  388. {
  389. b3Float4 dir;
  390. b3MprSupport_t v4;
  391. for (int i=0;i<B3_MPR_MAX_ITERATIONS;i++)
  392. //while (1)
  393. {
  394. // compute direction outside the portal (from v0 throught v1,v2,v3
  395. // face)
  396. b3PortalDir(portal, &dir);
  397. // test if origin is inside the portal
  398. if (portalEncapsulesOrigin(portal, &dir))
  399. return 0;
  400. // get next support point
  401. b3MprSupport(pairIndex,bodyIndexA,bodyIndexB,cpuBodyBuf,cpuConvexData,cpuCollidables,cpuVertices, sepAxis,&dir, &v4);
  402. // test if v4 can expand portal to contain origin and if portal
  403. // expanding doesn't reach given tolerance
  404. if (!portalCanEncapsuleOrigin(portal, &v4, &dir)
  405. || portalReachTolerance(portal, &v4, &dir))
  406. {
  407. return -1;
  408. }
  409. // v1-v2-v3 triangle must be rearranged to face outside Minkowski
  410. // difference (direction from v0).
  411. b3ExpandPortal(portal, &v4);
  412. }
  413. return -1;
  414. }
  415. B3_STATIC void b3FindPos(const b3MprSimplex_t *portal, b3Float4 *pos)
  416. {
  417. b3Float4 zero = b3MakeFloat4(0,0,0,0);
  418. b3Float4* b3mpr_vec3_origin = &zero;
  419. b3Float4 dir;
  420. size_t i;
  421. float b[4], sum, inv;
  422. b3Float4 vec, p1, p2;
  423. b3PortalDir(portal, &dir);
  424. // use barycentric coordinates of tetrahedron to find origin
  425. b3MprVec3Cross(&vec, &b3MprSimplexPoint(portal, 1)->v,
  426. &b3MprSimplexPoint(portal, 2)->v);
  427. b[0] = b3MprVec3Dot(&vec, &b3MprSimplexPoint(portal, 3)->v);
  428. b3MprVec3Cross(&vec, &b3MprSimplexPoint(portal, 3)->v,
  429. &b3MprSimplexPoint(portal, 2)->v);
  430. b[1] = b3MprVec3Dot(&vec, &b3MprSimplexPoint(portal, 0)->v);
  431. b3MprVec3Cross(&vec, &b3MprSimplexPoint(portal, 0)->v,
  432. &b3MprSimplexPoint(portal, 1)->v);
  433. b[2] = b3MprVec3Dot(&vec, &b3MprSimplexPoint(portal, 3)->v);
  434. b3MprVec3Cross(&vec, &b3MprSimplexPoint(portal, 2)->v,
  435. &b3MprSimplexPoint(portal, 1)->v);
  436. b[3] = b3MprVec3Dot(&vec, &b3MprSimplexPoint(portal, 0)->v);
  437. sum = b[0] + b[1] + b[2] + b[3];
  438. if (b3MprIsZero(sum) || sum < 0.f){
  439. b[0] = 0.f;
  440. b3MprVec3Cross(&vec, &b3MprSimplexPoint(portal, 2)->v,
  441. &b3MprSimplexPoint(portal, 3)->v);
  442. b[1] = b3MprVec3Dot(&vec, &dir);
  443. b3MprVec3Cross(&vec, &b3MprSimplexPoint(portal, 3)->v,
  444. &b3MprSimplexPoint(portal, 1)->v);
  445. b[2] = b3MprVec3Dot(&vec, &dir);
  446. b3MprVec3Cross(&vec, &b3MprSimplexPoint(portal, 1)->v,
  447. &b3MprSimplexPoint(portal, 2)->v);
  448. b[3] = b3MprVec3Dot(&vec, &dir);
  449. sum = b[1] + b[2] + b[3];
  450. }
  451. inv = 1.f / sum;
  452. b3MprVec3Copy(&p1, b3mpr_vec3_origin);
  453. b3MprVec3Copy(&p2, b3mpr_vec3_origin);
  454. for (i = 0; i < 4; i++){
  455. b3MprVec3Copy(&vec, &b3MprSimplexPoint(portal, i)->v1);
  456. b3MprVec3Scale(&vec, b[i]);
  457. b3MprVec3Add(&p1, &vec);
  458. b3MprVec3Copy(&vec, &b3MprSimplexPoint(portal, i)->v2);
  459. b3MprVec3Scale(&vec, b[i]);
  460. b3MprVec3Add(&p2, &vec);
  461. }
  462. b3MprVec3Scale(&p1, inv);
  463. b3MprVec3Scale(&p2, inv);
  464. b3MprVec3Copy(pos, &p1);
  465. b3MprVec3Add(pos, &p2);
  466. b3MprVec3Scale(pos, 0.5);
  467. }
  468. inline float b3MprVec3Dist2(const b3Float4 *a, const b3Float4 *b)
  469. {
  470. b3Float4 ab;
  471. b3MprVec3Sub2(&ab, a, b);
  472. return b3MprVec3Len2(&ab);
  473. }
  474. inline float _b3MprVec3PointSegmentDist2(const b3Float4 *P,
  475. const b3Float4 *x0,
  476. const b3Float4 *b,
  477. b3Float4 *witness)
  478. {
  479. // The computation comes from solving equation of segment:
  480. // S(t) = x0 + t.d
  481. // where - x0 is initial point of segment
  482. // - d is direction of segment from x0 (|d| > 0)
  483. // - t belongs to <0, 1> interval
  484. //
  485. // Than, distance from a segment to some point P can be expressed:
  486. // D(t) = |x0 + t.d - P|^2
  487. // which is distance from any point on segment. Minimization
  488. // of this function brings distance from P to segment.
  489. // Minimization of D(t) leads to simple quadratic equation that's
  490. // solving is straightforward.
  491. //
  492. // Bonus of this method is witness point for free.
  493. float dist, t;
  494. b3Float4 d, a;
  495. // direction of segment
  496. b3MprVec3Sub2(&d, b, x0);
  497. // precompute vector from P to x0
  498. b3MprVec3Sub2(&a, x0, P);
  499. t = -1.f * b3MprVec3Dot(&a, &d);
  500. t /= b3MprVec3Len2(&d);
  501. if (t < 0.f || b3MprIsZero(t)){
  502. dist = b3MprVec3Dist2(x0, P);
  503. if (witness)
  504. b3MprVec3Copy(witness, x0);
  505. }else if (t > 1.f || b3MprEq(t, 1.f)){
  506. dist = b3MprVec3Dist2(b, P);
  507. if (witness)
  508. b3MprVec3Copy(witness, b);
  509. }else{
  510. if (witness){
  511. b3MprVec3Copy(witness, &d);
  512. b3MprVec3Scale(witness, t);
  513. b3MprVec3Add(witness, x0);
  514. dist = b3MprVec3Dist2(witness, P);
  515. }else{
  516. // recycling variables
  517. b3MprVec3Scale(&d, t);
  518. b3MprVec3Add(&d, &a);
  519. dist = b3MprVec3Len2(&d);
  520. }
  521. }
  522. return dist;
  523. }
  524. inline float b3MprVec3PointTriDist2(const b3Float4 *P,
  525. const b3Float4 *x0, const b3Float4 *B,
  526. const b3Float4 *C,
  527. b3Float4 *witness)
  528. {
  529. // Computation comes from analytic expression for triangle (x0, B, C)
  530. // T(s, t) = x0 + s.d1 + t.d2, where d1 = B - x0 and d2 = C - x0 and
  531. // Then equation for distance is:
  532. // D(s, t) = | T(s, t) - P |^2
  533. // This leads to minimization of quadratic function of two variables.
  534. // The solution from is taken only if s is between 0 and 1, t is
  535. // between 0 and 1 and t + s < 1, otherwise distance from segment is
  536. // computed.
  537. b3Float4 d1, d2, a;
  538. float u, v, w, p, q, r;
  539. float s, t, dist, dist2;
  540. b3Float4 witness2;
  541. b3MprVec3Sub2(&d1, B, x0);
  542. b3MprVec3Sub2(&d2, C, x0);
  543. b3MprVec3Sub2(&a, x0, P);
  544. u = b3MprVec3Dot(&a, &a);
  545. v = b3MprVec3Dot(&d1, &d1);
  546. w = b3MprVec3Dot(&d2, &d2);
  547. p = b3MprVec3Dot(&a, &d1);
  548. q = b3MprVec3Dot(&a, &d2);
  549. r = b3MprVec3Dot(&d1, &d2);
  550. s = (q * r - w * p) / (w * v - r * r);
  551. t = (-s * r - q) / w;
  552. if ((b3MprIsZero(s) || s > 0.f)
  553. && (b3MprEq(s, 1.f) || s < 1.f)
  554. && (b3MprIsZero(t) || t > 0.f)
  555. && (b3MprEq(t, 1.f) || t < 1.f)
  556. && (b3MprEq(t + s, 1.f) || t + s < 1.f)){
  557. if (witness){
  558. b3MprVec3Scale(&d1, s);
  559. b3MprVec3Scale(&d2, t);
  560. b3MprVec3Copy(witness, x0);
  561. b3MprVec3Add(witness, &d1);
  562. b3MprVec3Add(witness, &d2);
  563. dist = b3MprVec3Dist2(witness, P);
  564. }else{
  565. dist = s * s * v;
  566. dist += t * t * w;
  567. dist += 2.f * s * t * r;
  568. dist += 2.f * s * p;
  569. dist += 2.f * t * q;
  570. dist += u;
  571. }
  572. }else{
  573. dist = _b3MprVec3PointSegmentDist2(P, x0, B, witness);
  574. dist2 = _b3MprVec3PointSegmentDist2(P, x0, C, &witness2);
  575. if (dist2 < dist){
  576. dist = dist2;
  577. if (witness)
  578. b3MprVec3Copy(witness, &witness2);
  579. }
  580. dist2 = _b3MprVec3PointSegmentDist2(P, B, C, &witness2);
  581. if (dist2 < dist){
  582. dist = dist2;
  583. if (witness)
  584. b3MprVec3Copy(witness, &witness2);
  585. }
  586. }
  587. return dist;
  588. }
  589. B3_STATIC void b3FindPenetr(int pairIndex,int bodyIndexA, int bodyIndexB, b3ConstArray(b3RigidBodyData_t) cpuBodyBuf,
  590. b3ConstArray(b3ConvexPolyhedronData_t) cpuConvexData,
  591. b3ConstArray(b3Collidable_t) cpuCollidables,
  592. b3ConstArray(b3Float4) cpuVertices,
  593. __global b3Float4* sepAxis,
  594. b3MprSimplex_t *portal,
  595. float *depth, b3Float4 *pdir, b3Float4 *pos)
  596. {
  597. b3Float4 dir;
  598. b3MprSupport_t v4;
  599. unsigned long iterations;
  600. b3Float4 zero = b3MakeFloat4(0,0,0,0);
  601. b3Float4* b3mpr_vec3_origin = &zero;
  602. iterations = 1UL;
  603. for (int i=0;i<B3_MPR_MAX_ITERATIONS;i++)
  604. //while (1)
  605. {
  606. // compute portal direction and obtain next support point
  607. b3PortalDir(portal, &dir);
  608. b3MprSupport(pairIndex,bodyIndexA,bodyIndexB,cpuBodyBuf,cpuConvexData,cpuCollidables,cpuVertices, sepAxis,&dir, &v4);
  609. // reached tolerance -> find penetration info
  610. if (portalReachTolerance(portal, &v4, &dir)
  611. || iterations ==B3_MPR_MAX_ITERATIONS)
  612. {
  613. *depth = b3MprVec3PointTriDist2(b3mpr_vec3_origin,&b3MprSimplexPoint(portal, 1)->v,&b3MprSimplexPoint(portal, 2)->v,&b3MprSimplexPoint(portal, 3)->v,pdir);
  614. *depth = B3_MPR_SQRT(*depth);
  615. if (b3MprIsZero((*pdir).x) && b3MprIsZero((*pdir).y) && b3MprIsZero((*pdir).z))
  616. {
  617. *pdir = dir;
  618. }
  619. b3MprVec3Normalize(pdir);
  620. // barycentric coordinates:
  621. b3FindPos(portal, pos);
  622. return;
  623. }
  624. b3ExpandPortal(portal, &v4);
  625. iterations++;
  626. }
  627. }
  628. B3_STATIC void b3FindPenetrTouch(b3MprSimplex_t *portal,float *depth, b3Float4 *dir, b3Float4 *pos)
  629. {
  630. // Touching contact on portal's v1 - so depth is zero and direction
  631. // is unimportant and pos can be guessed
  632. *depth = 0.f;
  633. b3Float4 zero = b3MakeFloat4(0,0,0,0);
  634. b3Float4* b3mpr_vec3_origin = &zero;
  635. b3MprVec3Copy(dir, b3mpr_vec3_origin);
  636. b3MprVec3Copy(pos, &b3MprSimplexPoint(portal, 1)->v1);
  637. b3MprVec3Add(pos, &b3MprSimplexPoint(portal, 1)->v2);
  638. b3MprVec3Scale(pos, 0.5);
  639. }
  640. B3_STATIC void b3FindPenetrSegment(b3MprSimplex_t *portal,
  641. float *depth, b3Float4 *dir, b3Float4 *pos)
  642. {
  643. // Origin lies on v0-v1 segment.
  644. // Depth is distance to v1, direction also and position must be
  645. // computed
  646. b3MprVec3Copy(pos, &b3MprSimplexPoint(portal, 1)->v1);
  647. b3MprVec3Add(pos, &b3MprSimplexPoint(portal, 1)->v2);
  648. b3MprVec3Scale(pos, 0.5f);
  649. b3MprVec3Copy(dir, &b3MprSimplexPoint(portal, 1)->v);
  650. *depth = B3_MPR_SQRT(b3MprVec3Len2(dir));
  651. b3MprVec3Normalize(dir);
  652. }
  653. inline int b3MprPenetration(int pairIndex, int bodyIndexA, int bodyIndexB,
  654. b3ConstArray(b3RigidBodyData_t) cpuBodyBuf,
  655. b3ConstArray(b3ConvexPolyhedronData_t) cpuConvexData,
  656. b3ConstArray(b3Collidable_t) cpuCollidables,
  657. b3ConstArray(b3Float4) cpuVertices,
  658. __global b3Float4* sepAxis,
  659. __global int* hasSepAxis,
  660. float *depthOut, b3Float4* dirOut, b3Float4* posOut)
  661. {
  662. b3MprSimplex_t portal;
  663. // if (!hasSepAxis[pairIndex])
  664. // return -1;
  665. hasSepAxis[pairIndex] = 0;
  666. int res;
  667. // Phase 1: Portal discovery
  668. res = b3DiscoverPortal(pairIndex,bodyIndexA,bodyIndexB,cpuBodyBuf,cpuConvexData,cpuCollidables,cpuVertices,sepAxis,hasSepAxis, &portal);
  669. //sepAxis[pairIndex] = *pdir;//or -dir?
  670. switch (res)
  671. {
  672. case 0:
  673. {
  674. // Phase 2: Portal refinement
  675. res = b3RefinePortal(pairIndex,bodyIndexA,bodyIndexB,cpuBodyBuf,cpuConvexData,cpuCollidables,cpuVertices, sepAxis,&portal);
  676. if (res < 0)
  677. return -1;
  678. // Phase 3. Penetration info
  679. b3FindPenetr(pairIndex,bodyIndexA,bodyIndexB,cpuBodyBuf,cpuConvexData,cpuCollidables,cpuVertices, sepAxis,&portal, depthOut, dirOut, posOut);
  680. hasSepAxis[pairIndex] = 1;
  681. sepAxis[pairIndex] = -*dirOut;
  682. break;
  683. }
  684. case 1:
  685. {
  686. // Touching contact on portal's v1.
  687. b3FindPenetrTouch(&portal, depthOut, dirOut, posOut);
  688. break;
  689. }
  690. case 2:
  691. {
  692. b3FindPenetrSegment( &portal, depthOut, dirOut, posOut);
  693. break;
  694. }
  695. default:
  696. {
  697. hasSepAxis[pairIndex]=0;
  698. //if (res < 0)
  699. //{
  700. // Origin isn't inside portal - no collision.
  701. return -1;
  702. //}
  703. }
  704. };
  705. return 0;
  706. };
  707. #endif //B3_MPR_PENETRATION_H