ConvexHullBuilder.cpp 41 KB

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