ConvexHullBuilder.cpp 41 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465
  1. // Jolt Physics Library (https://github.com/jrouwe/JoltPhysics)
  2. // SPDX-FileCopyrightText: 2021 Jorrit Rouwe
  3. // SPDX-License-Identifier: MIT
  4. #include <Jolt/Jolt.h>
  5. #include <Jolt/Geometry/ConvexHullBuilder.h>
  6. #include <Jolt/Geometry/ConvexHullBuilder2D.h>
  7. #include <Jolt/Geometry/ClosestPoint.h>
  8. #include <Jolt/Core/StringTools.h>
  9. #include <Jolt/Core/UnorderedSet.h>
  10. JPH_SUPPRESS_WARNINGS_STD_BEGIN
  11. #include <fstream>
  12. JPH_SUPPRESS_WARNINGS_STD_END
  13. #ifdef JPH_CONVEX_BUILDER_DEBUG
  14. #include <Jolt/Renderer/DebugRenderer.h>
  15. #endif
  16. JPH_NAMESPACE_BEGIN
  17. ConvexHullBuilder::Face::~Face()
  18. {
  19. // Free all edges
  20. Edge *e = mFirstEdge;
  21. if (e != nullptr)
  22. {
  23. do
  24. {
  25. Edge *next = e->mNextEdge;
  26. delete e;
  27. e = next;
  28. } while (e != mFirstEdge);
  29. }
  30. }
  31. void ConvexHullBuilder::Face::CalculateNormalAndCentroid(const Vec3 *inPositions)
  32. {
  33. // Get point that we use to construct a triangle fan
  34. Edge *e = mFirstEdge;
  35. Vec3 y0 = inPositions[e->mStartIdx];
  36. // Get the 2nd point
  37. e = e->mNextEdge;
  38. Vec3 y1 = inPositions[e->mStartIdx];
  39. // Start accumulating the centroid
  40. mCentroid = y0 + y1;
  41. int n = 2;
  42. // Start accumulating the normal
  43. mNormal = Vec3::sZero();
  44. // Loop over remaining edges accumulating normals in a triangle fan fashion
  45. for (e = e->mNextEdge; e != mFirstEdge; e = e->mNextEdge)
  46. {
  47. // Get the 3rd point
  48. Vec3 y2 = inPositions[e->mStartIdx];
  49. // Calculate edges (counter clockwise)
  50. Vec3 e0 = y1 - y0;
  51. Vec3 e1 = y2 - y1;
  52. Vec3 e2 = y0 - y2;
  53. // The best normal is calculated by using the two shortest edges
  54. // See: https://box2d.org/posts/2014/01/troublesome-triangle/
  55. // The difference in normals is most pronounced when one edge is much smaller than the others (in which case the others must have roughly the same length).
  56. // Therefore we can suffice by just picking the shortest from 2 edges and use that with the 3rd edge to calculate the normal.
  57. // We first check which of the edges is shorter: e1 or e2
  58. UVec4 e1_shorter_than_e2 = Vec4::sLess(e1.DotV4(e1), e2.DotV4(e2));
  59. // We calculate both normals and then select the one that had the shortest edge for our normal (this avoids branching)
  60. Vec3 normal_e01 = e0.Cross(e1);
  61. Vec3 normal_e02 = e2.Cross(e0);
  62. mNormal += Vec3::sSelect(normal_e02, normal_e01, e1_shorter_than_e2);
  63. // Accumulate centroid
  64. mCentroid += y2;
  65. n++;
  66. // Update y1 for next triangle
  67. y1 = y2;
  68. }
  69. // Finalize centroid
  70. mCentroid /= float(n);
  71. }
  72. void ConvexHullBuilder::Face::Initialize(int inIdx0, int inIdx1, int inIdx2, const Vec3 *inPositions)
  73. {
  74. JPH_ASSERT(mFirstEdge == nullptr);
  75. JPH_ASSERT(inIdx0 != inIdx1 && inIdx0 != inIdx2 && inIdx1 != inIdx2);
  76. // Create 3 edges
  77. Edge *e0 = new Edge(this, inIdx0);
  78. Edge *e1 = new Edge(this, inIdx1);
  79. Edge *e2 = new Edge(this, inIdx2);
  80. // Link edges
  81. e0->mNextEdge = e1;
  82. e1->mNextEdge = e2;
  83. e2->mNextEdge = e0;
  84. mFirstEdge = e0;
  85. CalculateNormalAndCentroid(inPositions);
  86. }
  87. ConvexHullBuilder::ConvexHullBuilder(const Positions &inPositions) :
  88. mPositions(inPositions)
  89. {
  90. #ifdef JPH_CONVEX_BUILDER_DEBUG
  91. mIteration = 0;
  92. // Center the drawing of the first hull around the origin and calculate the delta offset between states
  93. mOffset = RVec3::sZero();
  94. if (mPositions.empty())
  95. {
  96. // No hull will be generated
  97. mDelta = Vec3::sZero();
  98. }
  99. else
  100. {
  101. Vec3 maxv = Vec3::sReplicate(-FLT_MAX), minv = Vec3::sReplicate(FLT_MAX);
  102. for (Vec3 v : mPositions)
  103. {
  104. minv = Vec3::sMin(minv, v);
  105. maxv = Vec3::sMax(maxv, v);
  106. mOffset -= v;
  107. }
  108. mOffset /= Real(mPositions.size());
  109. mDelta = Vec3((maxv - minv).GetX() + 0.5f, 0, 0);
  110. mOffset += mDelta; // Don't start at origin, we're already drawing the final hull there
  111. }
  112. #endif
  113. }
  114. void ConvexHullBuilder::FreeFaces()
  115. {
  116. for (Face *f : mFaces)
  117. delete f;
  118. mFaces.clear();
  119. }
  120. void ConvexHullBuilder::GetFaceForPoint(Vec3Arg inPoint, const Faces &inFaces, Face *&outFace, float &outDistSq) const
  121. {
  122. outFace = nullptr;
  123. outDistSq = 0.0f;
  124. for (Face *f : inFaces)
  125. if (!f->mRemoved)
  126. {
  127. // Determine distance to face
  128. float dot = f->mNormal.Dot(inPoint - f->mCentroid);
  129. if (dot > 0.0f)
  130. {
  131. float dist_sq = dot * dot / f->mNormal.LengthSq();
  132. if (dist_sq > outDistSq)
  133. {
  134. outFace = f;
  135. outDistSq = dist_sq;
  136. }
  137. }
  138. }
  139. }
  140. float ConvexHullBuilder::GetDistanceToEdgeSq(Vec3Arg inPoint, const Face *inFace) const
  141. {
  142. bool all_inside = true;
  143. float edge_dist_sq = FLT_MAX;
  144. // Test if it is inside the edges of the polygon
  145. Edge *edge = inFace->mFirstEdge;
  146. Vec3 p1 = mPositions[edge->GetPreviousEdge()->mStartIdx];
  147. do
  148. {
  149. Vec3 p2 = mPositions[edge->mStartIdx];
  150. if ((p2 - p1).Cross(inPoint - p1).Dot(inFace->mNormal) < 0.0f)
  151. {
  152. // It is outside
  153. all_inside = false;
  154. // Measure distance to this edge
  155. uint32 s;
  156. edge_dist_sq = min(edge_dist_sq, ClosestPoint::GetClosestPointOnLine(p1 - inPoint, p2 - inPoint, s).LengthSq());
  157. }
  158. p1 = p2;
  159. edge = edge->mNextEdge;
  160. } while (edge != inFace->mFirstEdge);
  161. return all_inside? 0.0f : edge_dist_sq;
  162. }
  163. bool ConvexHullBuilder::AssignPointToFace(int inPositionIdx, const Faces &inFaces, float inToleranceSq)
  164. {
  165. Vec3 point = mPositions[inPositionIdx];
  166. // Find the face for which the point is furthest away
  167. Face *best_face;
  168. float best_dist_sq;
  169. GetFaceForPoint(point, inFaces, best_face, best_dist_sq);
  170. if (best_face != nullptr)
  171. {
  172. // Check if this point is within the tolerance margin to the plane
  173. if (best_dist_sq <= inToleranceSq)
  174. {
  175. // Check distance to edges
  176. float dist_to_edge_sq = GetDistanceToEdgeSq(point, best_face);
  177. if (dist_to_edge_sq > inToleranceSq)
  178. {
  179. // Point is outside of the face and too far away to discard
  180. mCoplanarList.push_back({ inPositionIdx, dist_to_edge_sq });
  181. }
  182. }
  183. else
  184. {
  185. // This point is in front of the face, add it to the conflict list
  186. if (best_dist_sq > best_face->mFurthestPointDistanceSq)
  187. {
  188. // This point is futher away than any others, update the distance and add point as last point
  189. best_face->mFurthestPointDistanceSq = best_dist_sq;
  190. best_face->mConflictList.push_back(inPositionIdx);
  191. }
  192. else
  193. {
  194. // Not the furthest point, add it as the before last point
  195. best_face->mConflictList.push_back(best_face->mConflictList.back());
  196. best_face->mConflictList[best_face->mConflictList.size() - 2] = inPositionIdx;
  197. }
  198. return true;
  199. }
  200. }
  201. return false;
  202. }
  203. float ConvexHullBuilder::DetermineCoplanarDistance() const
  204. {
  205. // Formula as per: Implementing Quickhull - Dirk Gregorius.
  206. Vec3 vmax = Vec3::sZero();
  207. for (Vec3 v : mPositions)
  208. vmax = Vec3::sMax(vmax, v.Abs());
  209. return 3.0f * FLT_EPSILON * (vmax.GetX() + vmax.GetY() + vmax.GetZ());
  210. }
  211. int ConvexHullBuilder::GetNumVerticesUsed() const
  212. {
  213. UnorderedSet<int> used_verts;
  214. for (Face *f : mFaces)
  215. {
  216. Edge *e = f->mFirstEdge;
  217. do
  218. {
  219. used_verts.insert(e->mStartIdx);
  220. e = e->mNextEdge;
  221. } while (e != f->mFirstEdge);
  222. }
  223. return (int)used_verts.size();
  224. }
  225. bool ConvexHullBuilder::ContainsFace(const Array<int> &inIndices) const
  226. {
  227. for (Face *f : mFaces)
  228. {
  229. Edge *e = f->mFirstEdge;
  230. Array<int>::const_iterator index = find(inIndices.begin(), inIndices.end(), e->mStartIdx);
  231. if (index != inIndices.end())
  232. {
  233. size_t matches = 0;
  234. do
  235. {
  236. // Check if index matches
  237. if (*index != e->mStartIdx)
  238. break;
  239. // Increment number of matches
  240. matches++;
  241. // Next index in list of inIndices
  242. index++;
  243. if (index == inIndices.end())
  244. index = inIndices.begin();
  245. // Next edge
  246. e = e->mNextEdge;
  247. } while (e != f->mFirstEdge);
  248. if (matches == inIndices.size())
  249. return true;
  250. }
  251. }
  252. return false;
  253. }
  254. ConvexHullBuilder::EResult ConvexHullBuilder::Initialize(int inMaxVertices, float inTolerance, const char *&outError)
  255. {
  256. // Free the faces possibly left over from an earlier hull
  257. FreeFaces();
  258. // Test that we have at least 3 points
  259. if (mPositions.size() < 3)
  260. {
  261. outError = "Need at least 3 points to make a hull";
  262. return EResult::TooFewPoints;
  263. }
  264. // Determine a suitable tolerance for detecting that points are coplanar
  265. float coplanar_tolerance_sq = Square(DetermineCoplanarDistance());
  266. // Increase desired tolerance if accuracy doesn't allow it
  267. float tolerance_sq = max(coplanar_tolerance_sq, Square(inTolerance));
  268. // Find point furthest from the origin
  269. int idx1 = -1;
  270. float max_dist_sq = -1.0f;
  271. for (int i = 0; i < (int)mPositions.size(); ++i)
  272. {
  273. float dist_sq = mPositions[i].LengthSq();
  274. if (dist_sq > max_dist_sq)
  275. {
  276. max_dist_sq = dist_sq;
  277. idx1 = i;
  278. }
  279. }
  280. JPH_ASSERT(idx1 >= 0);
  281. // Find point that is furthest away from this point
  282. int idx2 = -1;
  283. max_dist_sq = -1.0f;
  284. for (int i = 0; i < (int)mPositions.size(); ++i)
  285. if (i != idx1)
  286. {
  287. float dist_sq = (mPositions[i] - mPositions[idx1]).LengthSq();
  288. if (dist_sq > max_dist_sq)
  289. {
  290. max_dist_sq = dist_sq;
  291. idx2 = i;
  292. }
  293. }
  294. JPH_ASSERT(idx2 >= 0);
  295. // Find point that forms the biggest triangle
  296. int idx3 = -1;
  297. float best_triangle_area_sq = -1.0f;
  298. for (int i = 0; i < (int)mPositions.size(); ++i)
  299. if (i != idx1 && i != idx2)
  300. {
  301. float triangle_area_sq = (mPositions[idx1] - mPositions[i]).Cross(mPositions[idx2] - mPositions[i]).LengthSq();
  302. if (triangle_area_sq > best_triangle_area_sq)
  303. {
  304. best_triangle_area_sq = triangle_area_sq;
  305. idx3 = i;
  306. }
  307. }
  308. JPH_ASSERT(idx3 >= 0);
  309. if (best_triangle_area_sq < cMinTriangleAreaSq)
  310. {
  311. outError = "Could not find a suitable initial triangle because its area was too small";
  312. return EResult::Degenerate;
  313. }
  314. // Check if we have only 3 vertices
  315. if (mPositions.size() == 3)
  316. {
  317. // Create two triangles (back to back)
  318. Face *t1 = CreateTriangle(idx1, idx2, idx3);
  319. Face *t2 = CreateTriangle(idx1, idx3, idx2);
  320. // Link faces edges
  321. sLinkFace(t1->mFirstEdge, t2->mFirstEdge->mNextEdge->mNextEdge);
  322. sLinkFace(t1->mFirstEdge->mNextEdge, t2->mFirstEdge->mNextEdge);
  323. sLinkFace(t1->mFirstEdge->mNextEdge->mNextEdge, t2->mFirstEdge);
  324. #ifdef JPH_CONVEX_BUILDER_DEBUG
  325. // Draw current state
  326. DrawState();
  327. #endif
  328. return EResult::Success;
  329. }
  330. // Find point that forms the biggest tetrahedron
  331. Vec3 initial_plane_normal = (mPositions[idx2] - mPositions[idx1]).Cross(mPositions[idx3] - mPositions[idx1]).Normalized();
  332. Vec3 initial_plane_centroid = (mPositions[idx1] + mPositions[idx2] + mPositions[idx3]) / 3.0f;
  333. int idx4 = -1;
  334. float max_dist = 0.0f;
  335. for (int i = 0; i < (int)mPositions.size(); ++i)
  336. if (i != idx1 && i != idx2 && i != idx3)
  337. {
  338. float dist = (mPositions[i] - initial_plane_centroid).Dot(initial_plane_normal);
  339. if (abs(dist) > abs(max_dist))
  340. {
  341. max_dist = dist;
  342. idx4 = i;
  343. }
  344. }
  345. // Check if the hull is coplanar
  346. if (Square(max_dist) <= 25.0f * coplanar_tolerance_sq)
  347. {
  348. // First project all points in 2D space
  349. Vec3 base1 = initial_plane_normal.GetNormalizedPerpendicular();
  350. Vec3 base2 = initial_plane_normal.Cross(base1);
  351. Array<Vec3> positions_2d;
  352. positions_2d.reserve(mPositions.size());
  353. for (Vec3 v : mPositions)
  354. positions_2d.emplace_back(base1.Dot(v), base2.Dot(v), 0.0f);
  355. // Build hull
  356. Array<int> edges_2d;
  357. ConvexHullBuilder2D builder_2d(positions_2d);
  358. ConvexHullBuilder2D::EResult result = builder_2d.Initialize(idx1, idx2, idx3, inMaxVertices, inTolerance, edges_2d);
  359. // Create faces (back to back)
  360. Face *f1 = CreateFace();
  361. Face *f2 = CreateFace();
  362. // Create edges for face 1
  363. Array<Edge *> edges_f1;
  364. edges_f1.reserve(edges_2d.size());
  365. for (int start_idx : edges_2d)
  366. {
  367. Edge *edge = new Edge(f1, start_idx);
  368. if (edges_f1.empty())
  369. f1->mFirstEdge = edge;
  370. else
  371. edges_f1.back()->mNextEdge = edge;
  372. edges_f1.push_back(edge);
  373. }
  374. edges_f1.back()->mNextEdge = f1->mFirstEdge;
  375. // Create edges for face 2
  376. Array<Edge *> edges_f2;
  377. edges_f2.reserve(edges_2d.size());
  378. for (int i = (int)edges_2d.size() - 1; i >= 0; --i)
  379. {
  380. Edge *edge = new Edge(f2, edges_2d[i]);
  381. if (edges_f2.empty())
  382. f2->mFirstEdge = edge;
  383. else
  384. edges_f2.back()->mNextEdge = edge;
  385. edges_f2.push_back(edge);
  386. }
  387. edges_f2.back()->mNextEdge = f2->mFirstEdge;
  388. // Link edges
  389. for (size_t i = 0; i < edges_2d.size(); ++i)
  390. sLinkFace(edges_f1[i], edges_f2[(2 * edges_2d.size() - 2 - i) % edges_2d.size()]);
  391. // Calculate the plane for both faces
  392. f1->CalculateNormalAndCentroid(mPositions.data());
  393. f2->mNormal = -f1->mNormal;
  394. f2->mCentroid = f1->mCentroid;
  395. #ifdef JPH_CONVEX_BUILDER_DEBUG
  396. // Draw current state
  397. DrawState();
  398. #endif
  399. return result == ConvexHullBuilder2D::EResult::MaxVerticesReached? EResult::MaxVerticesReached : EResult::Success;
  400. }
  401. // Ensure the planes are facing outwards
  402. if (max_dist < 0.0f)
  403. swap(idx2, idx3);
  404. // Create tetrahedron
  405. Face *t1 = CreateTriangle(idx1, idx2, idx4);
  406. Face *t2 = CreateTriangle(idx2, idx3, idx4);
  407. Face *t3 = CreateTriangle(idx3, idx1, idx4);
  408. Face *t4 = CreateTriangle(idx1, idx3, idx2);
  409. // Link face edges
  410. sLinkFace(t1->mFirstEdge, t4->mFirstEdge->mNextEdge->mNextEdge);
  411. sLinkFace(t1->mFirstEdge->mNextEdge, t2->mFirstEdge->mNextEdge->mNextEdge);
  412. sLinkFace(t1->mFirstEdge->mNextEdge->mNextEdge, t3->mFirstEdge->mNextEdge);
  413. sLinkFace(t2->mFirstEdge, t4->mFirstEdge->mNextEdge);
  414. sLinkFace(t2->mFirstEdge->mNextEdge, t3->mFirstEdge->mNextEdge->mNextEdge);
  415. sLinkFace(t3->mFirstEdge, t4->mFirstEdge);
  416. // Build the initial conflict lists
  417. Faces faces { t1, t2, t3, t4 };
  418. for (int idx = 0; idx < (int)mPositions.size(); ++idx)
  419. if (idx != idx1 && idx != idx2 && idx != idx3 && idx != idx4)
  420. AssignPointToFace(idx, faces, tolerance_sq);
  421. #ifdef JPH_CONVEX_BUILDER_DEBUG
  422. // Draw current state including conflict list
  423. DrawState(true);
  424. // Increment iteration counter
  425. ++mIteration;
  426. #endif
  427. // Overestimate of the actual amount of vertices we use, for limiting the amount of vertices in the hull
  428. int num_vertices_used = 4;
  429. // Loop through the remainder of the points and add them
  430. for (;;)
  431. {
  432. // Find the face with the furthest point on it
  433. Face *face_with_furthest_point = nullptr;
  434. float furthest_dist_sq = 0.0f;
  435. for (Face *f : mFaces)
  436. if (f->mFurthestPointDistanceSq > furthest_dist_sq)
  437. {
  438. furthest_dist_sq = f->mFurthestPointDistanceSq;
  439. face_with_furthest_point = f;
  440. }
  441. int furthest_point_idx;
  442. if (face_with_furthest_point != nullptr)
  443. {
  444. // Take the furthest point
  445. furthest_point_idx = face_with_furthest_point->mConflictList.back();
  446. face_with_furthest_point->mConflictList.pop_back();
  447. }
  448. else if (!mCoplanarList.empty())
  449. {
  450. // Try to assign points to faces (this also recalculates the distance to the hull for the coplanar vertices)
  451. CoplanarList coplanar;
  452. mCoplanarList.swap(coplanar);
  453. bool added = false;
  454. for (const Coplanar &c : coplanar)
  455. added |= AssignPointToFace(c.mPositionIdx, mFaces, tolerance_sq);
  456. // If we were able to assign a point, loop again to pick it up
  457. if (added)
  458. continue;
  459. // If the coplanar list is empty, there are no points left and we're done
  460. if (mCoplanarList.empty())
  461. break;
  462. do
  463. {
  464. // Find the vertex that is furthest from the hull
  465. CoplanarList::size_type best_idx = 0;
  466. float best_dist_sq = mCoplanarList.front().mDistanceSq;
  467. for (CoplanarList::size_type idx = 1; idx < mCoplanarList.size(); ++idx)
  468. {
  469. const Coplanar &c = mCoplanarList[idx];
  470. if (c.mDistanceSq > best_dist_sq)
  471. {
  472. best_idx = idx;
  473. best_dist_sq = c.mDistanceSq;
  474. }
  475. }
  476. // Swap it to the end
  477. swap(mCoplanarList[best_idx], mCoplanarList.back());
  478. // Remove it
  479. furthest_point_idx = mCoplanarList.back().mPositionIdx;
  480. mCoplanarList.pop_back();
  481. // Find the face for which the point is furthest away
  482. GetFaceForPoint(mPositions[furthest_point_idx], mFaces, face_with_furthest_point, best_dist_sq);
  483. } while (!mCoplanarList.empty() && face_with_furthest_point == nullptr);
  484. if (face_with_furthest_point == nullptr)
  485. break;
  486. }
  487. else
  488. {
  489. // If there are no more vertices, we're done
  490. break;
  491. }
  492. // Check if we have a limit on the max vertices that we should produce
  493. if (num_vertices_used >= inMaxVertices)
  494. {
  495. // Count the actual amount of used vertices (we did not take the removal of any vertices into account)
  496. num_vertices_used = GetNumVerticesUsed();
  497. // Check if there are too many
  498. if (num_vertices_used >= inMaxVertices)
  499. return EResult::MaxVerticesReached;
  500. }
  501. // We're about to add another vertex
  502. ++num_vertices_used;
  503. // Add the point to the hull
  504. Faces new_faces;
  505. AddPoint(face_with_furthest_point, furthest_point_idx, coplanar_tolerance_sq, new_faces);
  506. // Redistribute points on conflict lists belonging to removed faces
  507. for (const Face *face : mFaces)
  508. if (face->mRemoved)
  509. for (int idx : face->mConflictList)
  510. AssignPointToFace(idx, new_faces, tolerance_sq);
  511. // Permanently delete faces that we removed in AddPoint()
  512. GarbageCollectFaces();
  513. #ifdef JPH_CONVEX_BUILDER_DEBUG
  514. // Draw state at the end of this step including conflict list
  515. DrawState(true);
  516. // Increment iteration counter
  517. ++mIteration;
  518. #endif
  519. }
  520. // Check if we are left with a hull. It is possible that hull building fails if the points are nearly coplanar.
  521. if (mFaces.size() < 2)
  522. {
  523. outError = "Too few faces in hull";
  524. return EResult::TooFewFaces;
  525. }
  526. return EResult::Success;
  527. }
  528. void ConvexHullBuilder::AddPoint(Face *inFacingFace, int inIdx, float inCoplanarToleranceSq, Faces &outNewFaces)
  529. {
  530. // Get position
  531. Vec3 pos = mPositions[inIdx];
  532. #ifdef JPH_CONVEX_BUILDER_DEBUG
  533. // Draw point to be added
  534. DebugRenderer::sInstance->DrawMarker(cDrawScale * (mOffset + pos), Color::sYellow, 0.1f);
  535. DebugRenderer::sInstance->DrawText3D(cDrawScale * (mOffset + pos), ConvertToString(inIdx), Color::sWhite);
  536. #endif
  537. #ifdef JPH_ENABLE_ASSERTS
  538. // Check if structure is intact
  539. ValidateFaces();
  540. #endif
  541. // Find edge of convex hull of faces that are not facing the new vertex
  542. FullEdges edges;
  543. FindEdge(inFacingFace, pos, edges);
  544. JPH_ASSERT(edges.size() >= 3);
  545. // Create new faces
  546. outNewFaces.reserve(edges.size());
  547. for (const FullEdge &e : edges)
  548. {
  549. JPH_ASSERT(e.mStartIdx != e.mEndIdx);
  550. Face *f = CreateTriangle(e.mStartIdx, e.mEndIdx, inIdx);
  551. outNewFaces.push_back(f);
  552. }
  553. // Link edges
  554. for (Faces::size_type i = 0; i < outNewFaces.size(); ++i)
  555. {
  556. sLinkFace(outNewFaces[i]->mFirstEdge, edges[i].mNeighbourEdge);
  557. sLinkFace(outNewFaces[i]->mFirstEdge->mNextEdge, outNewFaces[(i + 1) % outNewFaces.size()]->mFirstEdge->mNextEdge->mNextEdge);
  558. }
  559. // Loop on faces that were modified until nothing needs to be checked anymore
  560. Faces affected_faces = outNewFaces;
  561. while (!affected_faces.empty())
  562. {
  563. // Take the next face
  564. Face *face = affected_faces.back();
  565. affected_faces.pop_back();
  566. if (!face->mRemoved)
  567. {
  568. // Merge with neighbour if this is a degenerate face
  569. MergeDegenerateFace(face, affected_faces);
  570. // Merge with coplanar neighbours (or when the neighbour forms a concave edge)
  571. if (!face->mRemoved)
  572. MergeCoplanarOrConcaveFaces(face, inCoplanarToleranceSq, affected_faces);
  573. }
  574. }
  575. #ifdef JPH_ENABLE_ASSERTS
  576. // Check if structure is intact
  577. ValidateFaces();
  578. #endif
  579. }
  580. void ConvexHullBuilder::GarbageCollectFaces()
  581. {
  582. for (int i = (int)mFaces.size() - 1; i >= 0; --i)
  583. {
  584. Face *f = mFaces[i];
  585. if (f->mRemoved)
  586. {
  587. FreeFace(f);
  588. mFaces.erase(mFaces.begin() + i);
  589. }
  590. }
  591. }
  592. ConvexHullBuilder::Face *ConvexHullBuilder::CreateFace()
  593. {
  594. // Call provider to create face
  595. Face *f = new Face();
  596. #ifdef JPH_CONVEX_BUILDER_DEBUG
  597. // Remember iteration counter
  598. f->mIteration = mIteration;
  599. #endif
  600. // Add to list
  601. mFaces.push_back(f);
  602. return f;
  603. }
  604. ConvexHullBuilder::Face *ConvexHullBuilder::CreateTriangle(int inIdx1, int inIdx2, int inIdx3)
  605. {
  606. Face *f = CreateFace();
  607. f->Initialize(inIdx1, inIdx2, inIdx3, mPositions.data());
  608. return f;
  609. }
  610. void ConvexHullBuilder::FreeFace(Face *inFace)
  611. {
  612. JPH_ASSERT(inFace->mRemoved);
  613. #ifdef JPH_ENABLE_ASSERTS
  614. // Make sure that this face is not connected
  615. Edge *e = inFace->mFirstEdge;
  616. if (e != nullptr)
  617. do
  618. {
  619. JPH_ASSERT(e->mNeighbourEdge == nullptr);
  620. e = e->mNextEdge;
  621. } while (e != inFace->mFirstEdge);
  622. #endif
  623. // Free the face
  624. delete inFace;
  625. }
  626. void ConvexHullBuilder::sLinkFace(Edge *inEdge1, Edge *inEdge2)
  627. {
  628. // Check not connected yet
  629. JPH_ASSERT(inEdge1->mNeighbourEdge == nullptr);
  630. JPH_ASSERT(inEdge2->mNeighbourEdge == nullptr);
  631. JPH_ASSERT(inEdge1->mFace != inEdge2->mFace);
  632. // Check vertices match
  633. JPH_ASSERT(inEdge1->mStartIdx == inEdge2->mNextEdge->mStartIdx);
  634. JPH_ASSERT(inEdge2->mStartIdx == inEdge1->mNextEdge->mStartIdx);
  635. // Link up
  636. inEdge1->mNeighbourEdge = inEdge2;
  637. inEdge2->mNeighbourEdge = inEdge1;
  638. }
  639. void ConvexHullBuilder::sUnlinkFace(Face *inFace)
  640. {
  641. // Unlink from neighbours
  642. Edge *e = inFace->mFirstEdge;
  643. do
  644. {
  645. if (e->mNeighbourEdge != nullptr)
  646. {
  647. // Validate that neighbour points to us
  648. JPH_ASSERT(e->mNeighbourEdge->mNeighbourEdge == e);
  649. // Unlink
  650. e->mNeighbourEdge->mNeighbourEdge = nullptr;
  651. e->mNeighbourEdge = nullptr;
  652. }
  653. e = e->mNextEdge;
  654. } while (e != inFace->mFirstEdge);
  655. }
  656. void ConvexHullBuilder::FindEdge(Face *inFacingFace, Vec3Arg inVertex, FullEdges &outEdges) const
  657. {
  658. // Assert that we were given an empty array
  659. JPH_ASSERT(outEdges.empty());
  660. // Should start with a facing face
  661. JPH_ASSERT(inFacingFace->IsFacing(inVertex));
  662. // Flag as removed
  663. inFacingFace->mRemoved = true;
  664. // Instead of recursing, we build our own stack with the information we need
  665. struct StackEntry
  666. {
  667. Edge * mFirstEdge;
  668. Edge * mCurrentEdge;
  669. };
  670. constexpr int cMaxEdgeLength = 128;
  671. StackEntry stack[cMaxEdgeLength];
  672. int cur_stack_pos = 0;
  673. static_assert(alignof(Edge) >= 2, "Need lowest bit to indicate to tell if we completed the loop");
  674. // Start with the face / edge provided
  675. stack[0].mFirstEdge = inFacingFace->mFirstEdge;
  676. stack[0].mCurrentEdge = reinterpret_cast<Edge *>(reinterpret_cast<uintptr_t>(inFacingFace->mFirstEdge) | 1); // Set lowest bit of pointer to make it different from the first edge
  677. for (;;)
  678. {
  679. StackEntry &cur_entry = stack[cur_stack_pos];
  680. // Next edge
  681. Edge *raw_e = cur_entry.mCurrentEdge;
  682. Edge *e = reinterpret_cast<Edge *>(reinterpret_cast<uintptr_t>(raw_e) & ~uintptr_t(1)); // Remove the lowest bit which was used to indicate that this is the first edge we're testing
  683. cur_entry.mCurrentEdge = e->mNextEdge;
  684. // If we're back at the first edge we've completed the face and we're done
  685. if (raw_e == cur_entry.mFirstEdge)
  686. {
  687. // This face needs to be removed, unlink it now, caller will free
  688. sUnlinkFace(e->mFace);
  689. // Pop from stack
  690. if (--cur_stack_pos < 0)
  691. break;
  692. }
  693. else
  694. {
  695. // Visit neighbour face
  696. Edge *ne = e->mNeighbourEdge;
  697. if (ne != nullptr)
  698. {
  699. Face *n = ne->mFace;
  700. if (!n->mRemoved)
  701. {
  702. // Check if vertex is on the front side of this face
  703. if (n->IsFacing(inVertex))
  704. {
  705. // Vertex on front, this face needs to be removed
  706. n->mRemoved = true;
  707. // Add element to the stack of elements to visit
  708. cur_stack_pos++;
  709. JPH_ASSERT(cur_stack_pos < cMaxEdgeLength);
  710. StackEntry &new_entry = stack[cur_stack_pos];
  711. new_entry.mFirstEdge = ne;
  712. new_entry.mCurrentEdge = ne->mNextEdge; // We don't need to test this edge again since we came from it
  713. }
  714. else
  715. {
  716. // Vertex behind, keep edge
  717. FullEdge full;
  718. full.mNeighbourEdge = ne;
  719. full.mStartIdx = e->mStartIdx;
  720. full.mEndIdx = ne->mStartIdx;
  721. outEdges.push_back(full);
  722. }
  723. }
  724. }
  725. }
  726. }
  727. // Assert that we have a fully connected loop
  728. #ifdef JPH_ENABLE_ASSERTS
  729. for (int i = 0; i < (int)outEdges.size(); ++i)
  730. JPH_ASSERT(outEdges[i].mEndIdx == outEdges[(i + 1) % outEdges.size()].mStartIdx);
  731. #endif
  732. #ifdef JPH_CONVEX_BUILDER_DEBUG
  733. // Draw edge of facing faces
  734. for (int i = 0; i < (int)outEdges.size(); ++i)
  735. DebugRenderer::sInstance->DrawArrow(cDrawScale * (mOffset + mPositions[outEdges[i].mStartIdx]), cDrawScale * (mOffset + mPositions[outEdges[i].mEndIdx]), Color::sWhite, 0.01f);
  736. DrawState();
  737. #endif
  738. }
  739. void ConvexHullBuilder::MergeFaces(Edge *inEdge)
  740. {
  741. // Get the face
  742. Face *face = inEdge->mFace;
  743. // Find the previous and next edge
  744. Edge *next_edge = inEdge->mNextEdge;
  745. Edge *prev_edge = inEdge->GetPreviousEdge();
  746. // Get the other face
  747. Edge *other_edge = inEdge->mNeighbourEdge;
  748. Face *other_face = other_edge->mFace;
  749. // Check if attempting to merge with self
  750. JPH_ASSERT(face != other_face);
  751. #ifdef JPH_CONVEX_BUILDER_DEBUG
  752. DrawWireFace(face, Color::sGreen);
  753. DrawWireFace(other_face, Color::sRed);
  754. DrawState();
  755. #endif
  756. // Loop over the edges of the other face and make them belong to inFace
  757. Edge *edge = other_edge->mNextEdge;
  758. prev_edge->mNextEdge = edge;
  759. for (;;)
  760. {
  761. edge->mFace = face;
  762. if (edge->mNextEdge == other_edge)
  763. {
  764. // Terminate when we are back at other_edge
  765. edge->mNextEdge = next_edge;
  766. break;
  767. }
  768. edge = edge->mNextEdge;
  769. }
  770. // If the first edge happens to be inEdge we need to fix it because this edge is no longer part of the face.
  771. // Note that we replace it with the first edge of the merged face so that if the MergeFace function is called
  772. // from a loop that loops around the face that it will still terminate after visiting all edges once.
  773. if (face->mFirstEdge == inEdge)
  774. face->mFirstEdge = prev_edge->mNextEdge;
  775. // Free the edges
  776. delete inEdge;
  777. delete other_edge;
  778. // Mark the other face as removed
  779. other_face->mFirstEdge = nullptr;
  780. other_face->mRemoved = true;
  781. // Recalculate plane
  782. face->CalculateNormalAndCentroid(mPositions.data());
  783. // Merge conflict lists
  784. if (face->mFurthestPointDistanceSq > other_face->mFurthestPointDistanceSq)
  785. {
  786. // This face has a point that's further away, make sure it remains the last one as we add the other points to this faces list
  787. face->mConflictList.insert(face->mConflictList.end() - 1, other_face->mConflictList.begin(), other_face->mConflictList.end());
  788. }
  789. else
  790. {
  791. // The other face has a point that's furthest away, add that list at the end.
  792. face->mConflictList.insert(face->mConflictList.end(), other_face->mConflictList.begin(), other_face->mConflictList.end());
  793. face->mFurthestPointDistanceSq = other_face->mFurthestPointDistanceSq;
  794. }
  795. other_face->mConflictList.clear();
  796. #ifdef JPH_CONVEX_BUILDER_DEBUG
  797. DrawWireFace(face, Color::sWhite);
  798. DrawState();
  799. #endif
  800. }
  801. void ConvexHullBuilder::MergeDegenerateFace(Face *inFace, Faces &ioAffectedFaces)
  802. {
  803. // Check area of face
  804. if (inFace->mNormal.LengthSq() < cMinTriangleAreaSq)
  805. {
  806. // Find longest edge, since this face is a sliver this should keep the face convex
  807. float max_length_sq = 0.0f;
  808. Edge *longest_edge = nullptr;
  809. Edge *e = inFace->mFirstEdge;
  810. Vec3 p1 = mPositions[e->mStartIdx];
  811. do
  812. {
  813. Edge *next = e->mNextEdge;
  814. Vec3 p2 = mPositions[next->mStartIdx];
  815. float length_sq = (p2 - p1).LengthSq();
  816. if (length_sq >= max_length_sq)
  817. {
  818. max_length_sq = length_sq;
  819. longest_edge = e;
  820. }
  821. p1 = p2;
  822. e = next;
  823. } while (e != inFace->mFirstEdge);
  824. // Merge with face on longest edge
  825. MergeFaces(longest_edge);
  826. // Remove any invalid edges
  827. RemoveInvalidEdges(inFace, ioAffectedFaces);
  828. }
  829. }
  830. void ConvexHullBuilder::MergeCoplanarOrConcaveFaces(Face *inFace, float inCoplanarToleranceSq, Faces &ioAffectedFaces)
  831. {
  832. bool merged = false;
  833. Edge *edge = inFace->mFirstEdge;
  834. do
  835. {
  836. // Store next edge since this edge can be removed
  837. Edge *next_edge = edge->mNextEdge;
  838. // Test if centroid of one face is above plane of the other face by inCoplanarToleranceSq.
  839. // If so we need to merge other face into inFace.
  840. const Face *other_face = edge->mNeighbourEdge->mFace;
  841. Vec3 delta_centroid = other_face->mCentroid - inFace->mCentroid;
  842. float dist_other_face_centroid = inFace->mNormal.Dot(delta_centroid);
  843. float signed_dist_other_face_centroid_sq = abs(dist_other_face_centroid) * dist_other_face_centroid;
  844. float dist_face_centroid = -other_face->mNormal.Dot(delta_centroid);
  845. float signed_dist_face_centroid_sq = abs(dist_face_centroid) * dist_face_centroid;
  846. float face_normal_len_sq = inFace->mNormal.LengthSq();
  847. float other_face_normal_len_sq = other_face->mNormal.LengthSq();
  848. if ((signed_dist_other_face_centroid_sq > -inCoplanarToleranceSq * face_normal_len_sq
  849. || signed_dist_face_centroid_sq > -inCoplanarToleranceSq * other_face_normal_len_sq)
  850. && inFace->mNormal.Dot(other_face->mNormal) > 0.0f) // Never merge faces that are back to back
  851. {
  852. MergeFaces(edge);
  853. merged = true;
  854. }
  855. edge = next_edge;
  856. } while (edge != inFace->mFirstEdge);
  857. if (merged)
  858. RemoveInvalidEdges(inFace, ioAffectedFaces);
  859. }
  860. void ConvexHullBuilder::sMarkAffected(Face *inFace, Faces &ioAffectedFaces)
  861. {
  862. if (find(ioAffectedFaces.begin(), ioAffectedFaces.end(), inFace) == ioAffectedFaces.end())
  863. ioAffectedFaces.push_back(inFace);
  864. }
  865. void ConvexHullBuilder::RemoveInvalidEdges(Face *inFace, Faces &ioAffectedFaces)
  866. {
  867. // This marks that the plane needs to be recalculated (we delay this until the end of the
  868. // function since we don't use the plane and we want to avoid calculating it multiple times)
  869. bool recalculate_plane = false;
  870. // We keep going through this loop until no more edges were removed
  871. bool removed;
  872. do
  873. {
  874. removed = false;
  875. // Loop over all edges in this face
  876. Edge *edge = inFace->mFirstEdge;
  877. Face *neighbour_face = edge->mNeighbourEdge->mFace;
  878. do
  879. {
  880. Edge *next_edge = edge->mNextEdge;
  881. Face *next_neighbour_face = next_edge->mNeighbourEdge->mFace;
  882. if (neighbour_face == inFace)
  883. {
  884. // We only remove 1 edge at a time, check if this edge's next edge is our neighbour.
  885. // If this check fails, we will continue to scan along the edge until we find an edge where this is the case.
  886. if (edge->mNeighbourEdge == next_edge)
  887. {
  888. // This edge leads back to the starting point, this means the edge is interior and needs to be removed
  889. #ifdef JPH_CONVEX_BUILDER_DEBUG
  890. DrawWireFace(inFace, Color::sBlue);
  891. DrawState();
  892. #endif
  893. // Remove edge
  894. Edge *prev_edge = edge->GetPreviousEdge();
  895. prev_edge->mNextEdge = next_edge->mNextEdge;
  896. if (inFace->mFirstEdge == edge || inFace->mFirstEdge == next_edge)
  897. inFace->mFirstEdge = prev_edge;
  898. delete edge;
  899. delete next_edge;
  900. #ifdef JPH_CONVEX_BUILDER_DEBUG
  901. DrawWireFace(inFace, Color::sGreen);
  902. DrawState();
  903. #endif
  904. // Check if inFace now has only 2 edges left
  905. if (RemoveTwoEdgeFace(inFace, ioAffectedFaces))
  906. return; // Bail if face no longer exists
  907. // Restart the loop
  908. recalculate_plane = true;
  909. removed = true;
  910. break;
  911. }
  912. }
  913. else if (neighbour_face == next_neighbour_face)
  914. {
  915. // There are two edges that connect to the same face, we will remove the second one
  916. #ifdef JPH_CONVEX_BUILDER_DEBUG
  917. DrawWireFace(inFace, Color::sYellow);
  918. DrawWireFace(neighbour_face, Color::sRed);
  919. DrawState();
  920. #endif
  921. // First merge the neighbours edges
  922. Edge *neighbour_edge = next_edge->mNeighbourEdge;
  923. Edge *next_neighbour_edge = neighbour_edge->mNextEdge;
  924. if (neighbour_face->mFirstEdge == next_neighbour_edge)
  925. neighbour_face->mFirstEdge = neighbour_edge;
  926. neighbour_edge->mNextEdge = next_neighbour_edge->mNextEdge;
  927. neighbour_edge->mNeighbourEdge = edge;
  928. delete next_neighbour_edge;
  929. // Then merge my own edges
  930. if (inFace->mFirstEdge == next_edge)
  931. inFace->mFirstEdge = edge;
  932. edge->mNextEdge = next_edge->mNextEdge;
  933. edge->mNeighbourEdge = neighbour_edge;
  934. delete next_edge;
  935. #ifdef JPH_CONVEX_BUILDER_DEBUG
  936. DrawWireFace(inFace, Color::sYellow);
  937. DrawWireFace(neighbour_face, Color::sGreen);
  938. DrawState();
  939. #endif
  940. // Check if neighbour has only 2 edges left
  941. if (!RemoveTwoEdgeFace(neighbour_face, ioAffectedFaces))
  942. {
  943. // No, we need to recalculate its plane
  944. neighbour_face->CalculateNormalAndCentroid(mPositions.data());
  945. // Mark neighbour face as affected
  946. sMarkAffected(neighbour_face, ioAffectedFaces);
  947. }
  948. // Check if inFace now has only 2 edges left
  949. if (RemoveTwoEdgeFace(inFace, ioAffectedFaces))
  950. return; // Bail if face no longer exists
  951. // Restart loop
  952. recalculate_plane = true;
  953. removed = true;
  954. break;
  955. }
  956. // This edge is ok, go to the next edge
  957. edge = next_edge;
  958. neighbour_face = next_neighbour_face;
  959. } while (edge != inFace->mFirstEdge);
  960. } while (removed);
  961. // Recalculate plane?
  962. if (recalculate_plane)
  963. inFace->CalculateNormalAndCentroid(mPositions.data());
  964. }
  965. bool ConvexHullBuilder::RemoveTwoEdgeFace(Face *inFace, Faces &ioAffectedFaces) const
  966. {
  967. // Check if this face contains only 2 edges
  968. Edge *edge = inFace->mFirstEdge;
  969. Edge *next_edge = edge->mNextEdge;
  970. JPH_ASSERT(edge != next_edge); // 1 edge faces should not exist
  971. if (next_edge->mNextEdge == edge)
  972. {
  973. #ifdef JPH_CONVEX_BUILDER_DEBUG
  974. DrawWireFace(inFace, Color::sRed);
  975. DrawState();
  976. #endif
  977. // Schedule both neighbours for re-checking
  978. Edge *neighbour_edge = edge->mNeighbourEdge;
  979. Face *neighbour_face = neighbour_edge->mFace;
  980. Edge *next_neighbour_edge = next_edge->mNeighbourEdge;
  981. Face *next_neighbour_face = next_neighbour_edge->mFace;
  982. sMarkAffected(neighbour_face, ioAffectedFaces);
  983. sMarkAffected(next_neighbour_face, ioAffectedFaces);
  984. // Link my neighbours to each other
  985. neighbour_edge->mNeighbourEdge = next_neighbour_edge;
  986. next_neighbour_edge->mNeighbourEdge = neighbour_edge;
  987. // Unlink my edges
  988. edge->mNeighbourEdge = nullptr;
  989. next_edge->mNeighbourEdge = nullptr;
  990. // Mark this face as removed
  991. inFace->mRemoved = true;
  992. return true;
  993. }
  994. return false;
  995. }
  996. #ifdef JPH_ENABLE_ASSERTS
  997. void ConvexHullBuilder::DumpFace(const Face *inFace) const
  998. {
  999. Trace("f:0x%p", inFace);
  1000. const Edge *e = inFace->mFirstEdge;
  1001. do
  1002. {
  1003. Trace("\te:0x%p { i:%d e:0x%p f:0x%p }", e, e->mStartIdx, e->mNeighbourEdge, e->mNeighbourEdge->mFace);
  1004. e = e->mNextEdge;
  1005. } while (e != inFace->mFirstEdge);
  1006. }
  1007. void ConvexHullBuilder::DumpFaces() const
  1008. {
  1009. Trace("Dump Faces:");
  1010. for (const Face *f : mFaces)
  1011. if (!f->mRemoved)
  1012. DumpFace(f);
  1013. }
  1014. void ConvexHullBuilder::ValidateFace(const Face *inFace) const
  1015. {
  1016. if (inFace->mRemoved)
  1017. {
  1018. const Edge *e = inFace->mFirstEdge;
  1019. if (e != nullptr)
  1020. do
  1021. {
  1022. JPH_ASSERT(e->mNeighbourEdge == nullptr);
  1023. e = e->mNextEdge;
  1024. } while (e != inFace->mFirstEdge);
  1025. }
  1026. else
  1027. {
  1028. int edge_count = 0;
  1029. const Edge *e = inFace->mFirstEdge;
  1030. do
  1031. {
  1032. // Count edge
  1033. ++edge_count;
  1034. // Validate that adjacent faces are all different
  1035. if (mFaces.size() > 2)
  1036. for (const Edge *other_edge = e->mNextEdge; other_edge != inFace->mFirstEdge; other_edge = other_edge->mNextEdge)
  1037. JPH_ASSERT(e->mNeighbourEdge->mFace != other_edge->mNeighbourEdge->mFace);
  1038. // Assert that the face is correct
  1039. JPH_ASSERT(e->mFace == inFace);
  1040. // Assert that we have a neighbour
  1041. const Edge *nb_edge = e->mNeighbourEdge;
  1042. JPH_ASSERT(nb_edge != nullptr);
  1043. if (nb_edge != nullptr)
  1044. {
  1045. // Assert that our neighbours edge points to us
  1046. JPH_ASSERT(nb_edge->mNeighbourEdge == e);
  1047. // Assert that it belongs to a different face
  1048. JPH_ASSERT(nb_edge->mFace != inFace);
  1049. // Assert that the next edge of the neighbour points to the same vertex as this edge's vertex
  1050. JPH_ASSERT(nb_edge->mNextEdge->mStartIdx == e->mStartIdx);
  1051. // Assert that my next edge points to the same vertex as my neighbours vertex
  1052. JPH_ASSERT(e->mNextEdge->mStartIdx == nb_edge->mStartIdx);
  1053. }
  1054. e = e->mNextEdge;
  1055. } while (e != inFace->mFirstEdge);
  1056. // Assert that we have 3 or more edges
  1057. JPH_ASSERT(edge_count >= 3);
  1058. }
  1059. }
  1060. void ConvexHullBuilder::ValidateFaces() const
  1061. {
  1062. for (const Face *f : mFaces)
  1063. ValidateFace(f);
  1064. }
  1065. #endif // JPH_ENABLE_ASSERTS
  1066. void ConvexHullBuilder::GetCenterOfMassAndVolume(Vec3 &outCenterOfMass, float &outVolume) const
  1067. {
  1068. // Fourth point is the average of all face centroids
  1069. Vec3 v4 = Vec3::sZero();
  1070. for (const Face *f : mFaces)
  1071. v4 += f->mCentroid;
  1072. v4 /= float(mFaces.size());
  1073. // Calculate mass and center of mass of this convex hull by summing all tetrahedrons
  1074. outVolume = 0.0f;
  1075. outCenterOfMass = Vec3::sZero();
  1076. for (const Face *f : mFaces)
  1077. {
  1078. // Get the first vertex that we'll use to create a triangle fan
  1079. Edge *e = f->mFirstEdge;
  1080. Vec3 v1 = mPositions[e->mStartIdx];
  1081. // Get the second vertex
  1082. e = e->mNextEdge;
  1083. Vec3 v2 = mPositions[e->mStartIdx];
  1084. for (e = e->mNextEdge; e != f->mFirstEdge; e = e->mNextEdge)
  1085. {
  1086. // Fetch the last point of the triangle
  1087. Vec3 v3 = mPositions[e->mStartIdx];
  1088. // Calculate center of mass and mass of this tetrahedron,
  1089. // see: https://en.wikipedia.org/wiki/Tetrahedron#Volume
  1090. float volume_tetrahedron = (v1 - v4).Dot((v2 - v4).Cross(v3 - v4)); // Needs to be divided by 6, postpone this until the end of the loop
  1091. Vec3 center_of_mass_tetrahedron = v1 + v2 + v3 + v4; // Needs to be divided by 4, postpone this until the end of the loop
  1092. // Accumulate results
  1093. outVolume += volume_tetrahedron;
  1094. outCenterOfMass += volume_tetrahedron * center_of_mass_tetrahedron;
  1095. // Update v2 for next triangle
  1096. v2 = v3;
  1097. }
  1098. }
  1099. // Calculate center of mass, fall back to average point in case there is no volume (everything is on a plane in this case)
  1100. if (outVolume > FLT_EPSILON)
  1101. outCenterOfMass /= 4.0f * outVolume;
  1102. else
  1103. outCenterOfMass = v4;
  1104. outVolume /= 6.0f;
  1105. }
  1106. void ConvexHullBuilder::DetermineMaxError(Face *&outFaceWithMaxError, float &outMaxError, int &outMaxErrorPositionIdx, float &outCoplanarDistance) const
  1107. {
  1108. outCoplanarDistance = DetermineCoplanarDistance();
  1109. // This measures the distance from a polygon to the furthest point outside of the hull
  1110. float max_error = 0.0f;
  1111. Face *max_error_face = nullptr;
  1112. int max_error_point = -1;
  1113. for (int i = 0; i < (int)mPositions.size(); ++i)
  1114. {
  1115. Vec3 v = mPositions[i];
  1116. // This measures the closest edge from all faces to point v
  1117. // Note that we take the min of all faces since there may be multiple near coplanar faces so if we were to test this per face
  1118. // we may find that a point is outside of a polygon and mark it as an error, while it is actually inside a nearly coplanar
  1119. // polygon.
  1120. float min_edge_dist_sq = FLT_MAX;
  1121. Face *min_edge_dist_face = nullptr;
  1122. for (Face *f : mFaces)
  1123. {
  1124. // Check if point is on or in front of plane
  1125. float normal_len = f->mNormal.Length();
  1126. JPH_ASSERT(normal_len > 0.0f);
  1127. float plane_dist = f->mNormal.Dot(v - f->mCentroid) / normal_len;
  1128. if (plane_dist > -outCoplanarDistance)
  1129. {
  1130. // Check distance to the edges of this face
  1131. float edge_dist_sq = GetDistanceToEdgeSq(v, f);
  1132. if (edge_dist_sq < min_edge_dist_sq)
  1133. {
  1134. min_edge_dist_sq = edge_dist_sq;
  1135. min_edge_dist_face = f;
  1136. }
  1137. // If the point is inside the polygon and the point is in front of the plane, measure the distance
  1138. if (edge_dist_sq == 0.0f && plane_dist > max_error)
  1139. {
  1140. max_error = plane_dist;
  1141. max_error_face = f;
  1142. max_error_point = i;
  1143. }
  1144. }
  1145. }
  1146. // If the minimum distance to an edge is further than our current max error, we use that as max error
  1147. float min_edge_dist = sqrt(min_edge_dist_sq);
  1148. if (min_edge_dist_face != nullptr && min_edge_dist > max_error)
  1149. {
  1150. max_error = min_edge_dist;
  1151. max_error_face = min_edge_dist_face;
  1152. max_error_point = i;
  1153. }
  1154. }
  1155. outFaceWithMaxError = max_error_face;
  1156. outMaxError = max_error;
  1157. outMaxErrorPositionIdx = max_error_point;
  1158. }
  1159. #ifdef JPH_CONVEX_BUILDER_DEBUG
  1160. void ConvexHullBuilder::DrawState(bool inDrawConflictList) const
  1161. {
  1162. // Draw origin
  1163. DebugRenderer::sInstance->DrawMarker(cDrawScale * mOffset, Color::sRed, 0.2f);
  1164. int face_idx = 0;
  1165. // Draw faces
  1166. for (const Face *f : mFaces)
  1167. if (!f->mRemoved)
  1168. {
  1169. Color iteration_color = Color::sGetDistinctColor(f->mIteration);
  1170. Color face_color = Color::sGetDistinctColor(face_idx++);
  1171. // First point
  1172. const Edge *e = f->mFirstEdge;
  1173. RVec3 p1 = cDrawScale * (mOffset + mPositions[e->mStartIdx]);
  1174. // Second point
  1175. e = e->mNextEdge;
  1176. RVec3 p2 = cDrawScale * (mOffset + mPositions[e->mStartIdx]);
  1177. // First line
  1178. DebugRenderer::sInstance->DrawLine(p1, p2, Color::sGrey);
  1179. do
  1180. {
  1181. // Third point
  1182. e = e->mNextEdge;
  1183. RVec3 p3 = cDrawScale * (mOffset + mPositions[e->mStartIdx]);
  1184. DebugRenderer::sInstance->DrawTriangle(p1, p2, p3, iteration_color);
  1185. DebugRenderer::sInstance->DrawLine(p2, p3, Color::sGrey);
  1186. p2 = p3;
  1187. }
  1188. while (e != f->mFirstEdge);
  1189. // Draw normal
  1190. RVec3 centroid = cDrawScale * (mOffset + f->mCentroid);
  1191. DebugRenderer::sInstance->DrawArrow(centroid, centroid + f->mNormal.NormalizedOr(Vec3::sZero()), face_color, 0.01f);
  1192. // Draw conflict list
  1193. if (inDrawConflictList)
  1194. for (int idx : f->mConflictList)
  1195. DebugRenderer::sInstance->DrawMarker(cDrawScale * (mOffset + mPositions[idx]), face_color, 0.05f);
  1196. }
  1197. // Offset to the right
  1198. mOffset += mDelta;
  1199. }
  1200. void ConvexHullBuilder::DrawWireFace(const Face *inFace, ColorArg inColor) const
  1201. {
  1202. const Edge *e = inFace->mFirstEdge;
  1203. RVec3 prev = cDrawScale * (mOffset + mPositions[e->mStartIdx]);
  1204. do
  1205. {
  1206. const Edge *next = e->mNextEdge;
  1207. RVec3 cur = cDrawScale * (mOffset + mPositions[next->mStartIdx]);
  1208. DebugRenderer::sInstance->DrawArrow(prev, cur, inColor, 0.01f);
  1209. DebugRenderer::sInstance->DrawText3D(prev, ConvertToString(e->mStartIdx), inColor);
  1210. e = next;
  1211. prev = cur;
  1212. } while (e != inFace->mFirstEdge);
  1213. }
  1214. void ConvexHullBuilder::DrawEdge(const Edge *inEdge, ColorArg inColor) const
  1215. {
  1216. RVec3 p1 = cDrawScale * (mOffset + mPositions[inEdge->mStartIdx]);
  1217. RVec3 p2 = cDrawScale * (mOffset + mPositions[inEdge->mNextEdge->mStartIdx]);
  1218. DebugRenderer::sInstance->DrawArrow(p1, p2, inColor, 0.01f);
  1219. }
  1220. #endif // JPH_CONVEX_BUILDER_DEBUG
  1221. #ifdef JPH_CONVEX_BUILDER_DUMP_SHAPE
  1222. void ConvexHullBuilder::DumpShape() const
  1223. {
  1224. static atomic<int> sShapeNo = 1;
  1225. int shape_no = sShapeNo++;
  1226. std::ofstream f;
  1227. f.open(StringFormat("dumped_shape%d.cpp", shape_no).c_str(), std::ofstream::out | std::ofstream::trunc);
  1228. if (!f.is_open())
  1229. return;
  1230. f << "{\n";
  1231. for (Vec3 v : mPositions)
  1232. f << StringFormat("\tVec3(%.9gf, %.9gf, %.9gf),\n", (double)v.GetX(), (double)v.GetY(), (double)v.GetZ());
  1233. f << "},\n";
  1234. }
  1235. #endif // JPH_CONVEX_BUILDER_DUMP_SHAPE
  1236. JPH_NAMESPACE_END