ConvexHullBuilder.cpp 41 KB

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