ConvexHullBuilder.cpp 41 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298129913001301130213031304130513061307130813091310131113121313131413151316131713181319132013211322132313241325132613271328132913301331133213331334133513361337133813391340134113421343134413451346134713481349135013511352135313541355135613571358135913601361136213631364136513661367136813691370137113721373137413751376137713781379138013811382138313841385138613871388138913901391139213931394139513961397139813991400140114021403140414051406140714081409141014111412141314141415141614171418141914201421142214231424142514261427142814291430143114321433143414351436143714381439144014411442144314441445144614471448144914501451145214531454145514561457
  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. return EResult::Success;
  520. }
  521. void ConvexHullBuilder::AddPoint(Face *inFacingFace, int inIdx, float inCoplanarToleranceSq, Faces &outNewFaces)
  522. {
  523. // Get position
  524. Vec3 pos = mPositions[inIdx];
  525. #ifdef JPH_CONVEX_BUILDER_DEBUG
  526. // Draw point to be added
  527. DebugRenderer::sInstance->DrawMarker(cDrawScale * (mOffset + pos), Color::sYellow, 0.1f);
  528. DebugRenderer::sInstance->DrawText3D(cDrawScale * (mOffset + pos), ConvertToString(inIdx), Color::sWhite);
  529. #endif
  530. #ifdef JPH_ENABLE_ASSERTS
  531. // Check if structure is intact
  532. ValidateFaces();
  533. #endif
  534. // Find edge of convex hull of faces that are not facing the new vertex
  535. FullEdges edges;
  536. FindEdge(inFacingFace, pos, edges);
  537. JPH_ASSERT(edges.size() >= 3);
  538. // Create new faces
  539. outNewFaces.reserve(edges.size());
  540. for (const FullEdge &e : edges)
  541. {
  542. JPH_ASSERT(e.mStartIdx != e.mEndIdx);
  543. Face *f = CreateTriangle(e.mStartIdx, e.mEndIdx, inIdx);
  544. outNewFaces.push_back(f);
  545. }
  546. // Link edges
  547. for (Faces::size_type i = 0; i < outNewFaces.size(); ++i)
  548. {
  549. sLinkFace(outNewFaces[i]->mFirstEdge, edges[i].mNeighbourEdge);
  550. sLinkFace(outNewFaces[i]->mFirstEdge->mNextEdge, outNewFaces[(i + 1) % outNewFaces.size()]->mFirstEdge->mNextEdge->mNextEdge);
  551. }
  552. // Loop on faces that were modified until nothing needs to be checked anymore
  553. Faces affected_faces = outNewFaces;
  554. while (!affected_faces.empty())
  555. {
  556. // Take the next face
  557. Face *face = affected_faces.back();
  558. affected_faces.pop_back();
  559. if (!face->mRemoved)
  560. {
  561. // Merge with neighbour if this is a degenerate face
  562. MergeDegenerateFace(face, affected_faces);
  563. // Merge with coplanar neighbours (or when the neighbour forms a concave edge)
  564. if (!face->mRemoved)
  565. MergeCoplanarOrConcaveFaces(face, inCoplanarToleranceSq, affected_faces);
  566. }
  567. }
  568. #ifdef JPH_ENABLE_ASSERTS
  569. // Check if structure is intact
  570. ValidateFaces();
  571. #endif
  572. }
  573. void ConvexHullBuilder::GarbageCollectFaces()
  574. {
  575. for (int i = (int)mFaces.size() - 1; i >= 0; --i)
  576. {
  577. Face *f = mFaces[i];
  578. if (f->mRemoved)
  579. {
  580. FreeFace(f);
  581. mFaces.erase(mFaces.begin() + i);
  582. }
  583. }
  584. }
  585. ConvexHullBuilder::Face *ConvexHullBuilder::CreateFace()
  586. {
  587. // Call provider to create face
  588. Face *f = new Face();
  589. #ifdef JPH_CONVEX_BUILDER_DEBUG
  590. // Remember iteration counter
  591. f->mIteration = mIteration;
  592. #endif
  593. // Add to list
  594. mFaces.push_back(f);
  595. return f;
  596. }
  597. ConvexHullBuilder::Face *ConvexHullBuilder::CreateTriangle(int inIdx1, int inIdx2, int inIdx3)
  598. {
  599. Face *f = CreateFace();
  600. f->Initialize(inIdx1, inIdx2, inIdx3, mPositions.data());
  601. return f;
  602. }
  603. void ConvexHullBuilder::FreeFace(Face *inFace)
  604. {
  605. JPH_ASSERT(inFace->mRemoved);
  606. #ifdef JPH_ENABLE_ASSERTS
  607. // Make sure that this face is not connected
  608. Edge *e = inFace->mFirstEdge;
  609. if (e != nullptr)
  610. do
  611. {
  612. JPH_ASSERT(e->mNeighbourEdge == nullptr);
  613. e = e->mNextEdge;
  614. } while (e != inFace->mFirstEdge);
  615. #endif
  616. // Free the face
  617. delete inFace;
  618. }
  619. void ConvexHullBuilder::sLinkFace(Edge *inEdge1, Edge *inEdge2)
  620. {
  621. // Check not connected yet
  622. JPH_ASSERT(inEdge1->mNeighbourEdge == nullptr);
  623. JPH_ASSERT(inEdge2->mNeighbourEdge == nullptr);
  624. JPH_ASSERT(inEdge1->mFace != inEdge2->mFace);
  625. // Check vertices match
  626. JPH_ASSERT(inEdge1->mStartIdx == inEdge2->mNextEdge->mStartIdx);
  627. JPH_ASSERT(inEdge2->mStartIdx == inEdge1->mNextEdge->mStartIdx);
  628. // Link up
  629. inEdge1->mNeighbourEdge = inEdge2;
  630. inEdge2->mNeighbourEdge = inEdge1;
  631. }
  632. void ConvexHullBuilder::sUnlinkFace(Face *inFace)
  633. {
  634. // Unlink from neighbours
  635. Edge *e = inFace->mFirstEdge;
  636. do
  637. {
  638. if (e->mNeighbourEdge != nullptr)
  639. {
  640. // Validate that neighbour points to us
  641. JPH_ASSERT(e->mNeighbourEdge->mNeighbourEdge == e);
  642. // Unlink
  643. e->mNeighbourEdge->mNeighbourEdge = nullptr;
  644. e->mNeighbourEdge = nullptr;
  645. }
  646. e = e->mNextEdge;
  647. } while (e != inFace->mFirstEdge);
  648. }
  649. void ConvexHullBuilder::FindEdge(Face *inFacingFace, Vec3Arg inVertex, FullEdges &outEdges) const
  650. {
  651. // Assert that we were given an empty array
  652. JPH_ASSERT(outEdges.empty());
  653. // Should start with a facing face
  654. JPH_ASSERT(inFacingFace->IsFacing(inVertex));
  655. // Flag as removed
  656. inFacingFace->mRemoved = true;
  657. // Instead of recursing, we build our own stack with the information we need
  658. struct StackEntry
  659. {
  660. Edge * mFirstEdge;
  661. Edge * mCurrentEdge;
  662. };
  663. constexpr int cMaxEdgeLength = 128;
  664. StackEntry stack[cMaxEdgeLength];
  665. int cur_stack_pos = 0;
  666. static_assert(alignof(Edge) >= 2, "Need lowest bit to indicate to tell if we completed the loop");
  667. // Start with the face / edge provided
  668. stack[0].mFirstEdge = inFacingFace->mFirstEdge;
  669. 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
  670. for (;;)
  671. {
  672. StackEntry &cur_entry = stack[cur_stack_pos];
  673. // Next edge
  674. Edge *raw_e = cur_entry.mCurrentEdge;
  675. 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
  676. cur_entry.mCurrentEdge = e->mNextEdge;
  677. // If we're back at the first edge we've completed the face and we're done
  678. if (raw_e == cur_entry.mFirstEdge)
  679. {
  680. // This face needs to be removed, unlink it now, caller will free
  681. sUnlinkFace(e->mFace);
  682. // Pop from stack
  683. if (--cur_stack_pos < 0)
  684. break;
  685. }
  686. else
  687. {
  688. // Visit neighbour face
  689. Edge *ne = e->mNeighbourEdge;
  690. if (ne != nullptr)
  691. {
  692. Face *n = ne->mFace;
  693. if (!n->mRemoved)
  694. {
  695. // Check if vertex is on the front side of this face
  696. if (n->IsFacing(inVertex))
  697. {
  698. // Vertex on front, this face needs to be removed
  699. n->mRemoved = true;
  700. // Add element to the stack of elements to visit
  701. cur_stack_pos++;
  702. JPH_ASSERT(cur_stack_pos < cMaxEdgeLength);
  703. StackEntry &new_entry = stack[cur_stack_pos];
  704. new_entry.mFirstEdge = ne;
  705. new_entry.mCurrentEdge = ne->mNextEdge; // We don't need to test this edge again since we came from it
  706. }
  707. else
  708. {
  709. // Vertex behind, keep edge
  710. FullEdge full;
  711. full.mNeighbourEdge = ne;
  712. full.mStartIdx = e->mStartIdx;
  713. full.mEndIdx = ne->mStartIdx;
  714. outEdges.push_back(full);
  715. }
  716. }
  717. }
  718. }
  719. }
  720. // Assert that we have a fully connected loop
  721. #ifdef JPH_ENABLE_ASSERTS
  722. for (int i = 0; i < (int)outEdges.size(); ++i)
  723. JPH_ASSERT(outEdges[i].mEndIdx == outEdges[(i + 1) % outEdges.size()].mStartIdx);
  724. #endif
  725. #ifdef JPH_CONVEX_BUILDER_DEBUG
  726. // Draw edge of facing faces
  727. for (int i = 0; i < (int)outEdges.size(); ++i)
  728. DebugRenderer::sInstance->DrawArrow(cDrawScale * (mPositions[outEdges[i].mStartIdx] + mOffset), cDrawScale * (mPositions[outEdges[i].mEndIdx] + mOffset), Color::sWhite, 0.01f);
  729. DrawState();
  730. #endif
  731. }
  732. void ConvexHullBuilder::MergeFaces(Edge *inEdge)
  733. {
  734. // Get the face
  735. Face *face = inEdge->mFace;
  736. // Find the previous and next edge
  737. Edge *next_edge = inEdge->mNextEdge;
  738. Edge *prev_edge = inEdge->GetPreviousEdge();
  739. // Get the other face
  740. Edge *other_edge = inEdge->mNeighbourEdge;
  741. Face *other_face = other_edge->mFace;
  742. // Check if attempting to merge with self
  743. JPH_ASSERT(face != other_face);
  744. #ifdef JPH_CONVEX_BUILDER_DEBUG
  745. DrawWireFace(face, Color::sGreen);
  746. DrawWireFace(other_face, Color::sRed);
  747. DrawState();
  748. #endif
  749. // Loop over the edges of the other face and make them belong to inFace
  750. Edge *edge = other_edge->mNextEdge;
  751. prev_edge->mNextEdge = edge;
  752. for (;;)
  753. {
  754. edge->mFace = face;
  755. if (edge->mNextEdge == other_edge)
  756. {
  757. // Terminate when we are back at other_edge
  758. edge->mNextEdge = next_edge;
  759. break;
  760. }
  761. edge = edge->mNextEdge;
  762. }
  763. // If the first edge happens to be inEdge we need to fix it because this edge is no longer part of the face.
  764. // Note that we replace it with the first edge of the merged face so that if the MergeFace function is called
  765. // from a loop that loops around the face that it will still terminate after visiting all edges once.
  766. if (face->mFirstEdge == inEdge)
  767. face->mFirstEdge = prev_edge->mNextEdge;
  768. // Free the edges
  769. delete inEdge;
  770. delete other_edge;
  771. // Mark the other face as removed
  772. other_face->mFirstEdge = nullptr;
  773. other_face->mRemoved = true;
  774. // Recalculate plane
  775. face->CalculateNormalAndCentroid(mPositions.data());
  776. // Merge conflict lists
  777. if (face->mFurthestPointDistanceSq > other_face->mFurthestPointDistanceSq)
  778. {
  779. // 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
  780. face->mConflictList.insert(face->mConflictList.end() - 1, other_face->mConflictList.begin(), other_face->mConflictList.end());
  781. }
  782. else
  783. {
  784. // The other face has a point that's furthest away, add that list at the end.
  785. face->mConflictList.insert(face->mConflictList.end(), other_face->mConflictList.begin(), other_face->mConflictList.end());
  786. face->mFurthestPointDistanceSq = other_face->mFurthestPointDistanceSq;
  787. }
  788. other_face->mConflictList.clear();
  789. #ifdef JPH_CONVEX_BUILDER_DEBUG
  790. DrawWireFace(face, Color::sWhite);
  791. DrawState();
  792. #endif
  793. }
  794. void ConvexHullBuilder::MergeDegenerateFace(Face *inFace, Faces &ioAffectedFaces)
  795. {
  796. // Check area of face
  797. if (inFace->mNormal.LengthSq() < cMinTriangleAreaSq)
  798. {
  799. // Find longest edge, since this face is a sliver this should keep the face convex
  800. float max_length_sq = 0.0f;
  801. Edge *longest_edge = nullptr;
  802. Edge *e = inFace->mFirstEdge;
  803. Vec3 p1 = mPositions[e->mStartIdx];
  804. do
  805. {
  806. Edge *next = e->mNextEdge;
  807. Vec3 p2 = mPositions[next->mStartIdx];
  808. float length_sq = (p2 - p1).LengthSq();
  809. if (length_sq >= max_length_sq)
  810. {
  811. max_length_sq = length_sq;
  812. longest_edge = e;
  813. }
  814. p1 = p2;
  815. e = next;
  816. } while (e != inFace->mFirstEdge);
  817. // Merge with face on longest edge
  818. MergeFaces(longest_edge);
  819. // Remove any invalid edges
  820. RemoveInvalidEdges(inFace, ioAffectedFaces);
  821. }
  822. }
  823. void ConvexHullBuilder::MergeCoplanarOrConcaveFaces(Face *inFace, float inCoplanarToleranceSq, Faces &ioAffectedFaces)
  824. {
  825. bool merged = false;
  826. Edge *edge = inFace->mFirstEdge;
  827. do
  828. {
  829. // Store next edge since this edge can be removed
  830. Edge *next_edge = edge->mNextEdge;
  831. // Test if centroid of one face is above plane of the other face by inCoplanarToleranceSq.
  832. // If so we need to merge other face into inFace.
  833. const Face *other_face = edge->mNeighbourEdge->mFace;
  834. Vec3 delta_centroid = other_face->mCentroid - inFace->mCentroid;
  835. float dist_other_face_centroid = inFace->mNormal.Dot(delta_centroid);
  836. float signed_dist_other_face_centroid_sq = abs(dist_other_face_centroid) * dist_other_face_centroid;
  837. float dist_face_centroid = -other_face->mNormal.Dot(delta_centroid);
  838. float signed_dist_face_centroid_sq = abs(dist_face_centroid) * dist_face_centroid;
  839. float face_normal_len_sq = inFace->mNormal.LengthSq();
  840. float other_face_normal_len_sq = other_face->mNormal.LengthSq();
  841. if ((signed_dist_other_face_centroid_sq > -inCoplanarToleranceSq * face_normal_len_sq
  842. || signed_dist_face_centroid_sq > -inCoplanarToleranceSq * other_face_normal_len_sq)
  843. && inFace->mNormal.Dot(other_face->mNormal) > 0.0f) // Never merge faces that are back to back
  844. {
  845. MergeFaces(edge);
  846. merged = true;
  847. }
  848. edge = next_edge;
  849. } while (edge != inFace->mFirstEdge);
  850. if (merged)
  851. RemoveInvalidEdges(inFace, ioAffectedFaces);
  852. }
  853. void ConvexHullBuilder::sMarkAffected(Face *inFace, Faces &ioAffectedFaces)
  854. {
  855. if (find(ioAffectedFaces.begin(), ioAffectedFaces.end(), inFace) == ioAffectedFaces.end())
  856. ioAffectedFaces.push_back(inFace);
  857. }
  858. void ConvexHullBuilder::RemoveInvalidEdges(Face *inFace, Faces &ioAffectedFaces)
  859. {
  860. // This marks that the plane needs to be recalculated (we delay this until the end of the
  861. // function since we don't use the plane and we want to avoid calculating it multiple times)
  862. bool recalculate_plane = false;
  863. // We keep going through this loop until no more edges were removed
  864. bool removed;
  865. do
  866. {
  867. removed = false;
  868. // Loop over all edges in this face
  869. Edge *edge = inFace->mFirstEdge;
  870. Face *neighbour_face = edge->mNeighbourEdge->mFace;
  871. do
  872. {
  873. Edge *next_edge = edge->mNextEdge;
  874. Face *next_neighbour_face = next_edge->mNeighbourEdge->mFace;
  875. if (neighbour_face == inFace)
  876. {
  877. // We only remove 1 edge at a time, check if this edge's next edge is our neighbour.
  878. // If this check fails, we will continue to scan along the edge until we find an edge where this is the case.
  879. if (edge->mNeighbourEdge == next_edge)
  880. {
  881. // This edge leads back to the starting point, this means the edge is interior and needs to be removed
  882. #ifdef JPH_CONVEX_BUILDER_DEBUG
  883. DrawWireFace(inFace, Color::sBlue);
  884. DrawState();
  885. #endif
  886. // Remove edge
  887. Edge *prev_edge = edge->GetPreviousEdge();
  888. prev_edge->mNextEdge = next_edge->mNextEdge;
  889. if (inFace->mFirstEdge == edge || inFace->mFirstEdge == next_edge)
  890. inFace->mFirstEdge = prev_edge;
  891. delete edge;
  892. delete next_edge;
  893. #ifdef JPH_CONVEX_BUILDER_DEBUG
  894. DrawWireFace(inFace, Color::sGreen);
  895. DrawState();
  896. #endif
  897. // Check if inFace now has only 2 edges left
  898. if (RemoveTwoEdgeFace(inFace, ioAffectedFaces))
  899. return; // Bail if face no longer exists
  900. // Restart the loop
  901. recalculate_plane = true;
  902. removed = true;
  903. break;
  904. }
  905. }
  906. else if (neighbour_face == next_neighbour_face)
  907. {
  908. // There are two edges that connect to the same face, we will remove the second one
  909. #ifdef JPH_CONVEX_BUILDER_DEBUG
  910. DrawWireFace(inFace, Color::sYellow);
  911. DrawWireFace(neighbour_face, Color::sRed);
  912. DrawState();
  913. #endif
  914. // First merge the neighbours edges
  915. Edge *neighbour_edge = next_edge->mNeighbourEdge;
  916. Edge *next_neighbour_edge = neighbour_edge->mNextEdge;
  917. if (neighbour_face->mFirstEdge == next_neighbour_edge)
  918. neighbour_face->mFirstEdge = neighbour_edge;
  919. neighbour_edge->mNextEdge = next_neighbour_edge->mNextEdge;
  920. neighbour_edge->mNeighbourEdge = edge;
  921. delete next_neighbour_edge;
  922. // Then merge my own edges
  923. if (inFace->mFirstEdge == next_edge)
  924. inFace->mFirstEdge = edge;
  925. edge->mNextEdge = next_edge->mNextEdge;
  926. edge->mNeighbourEdge = neighbour_edge;
  927. delete next_edge;
  928. #ifdef JPH_CONVEX_BUILDER_DEBUG
  929. DrawWireFace(inFace, Color::sYellow);
  930. DrawWireFace(neighbour_face, Color::sGreen);
  931. DrawState();
  932. #endif
  933. // Check if neighbour has only 2 edges left
  934. if (!RemoveTwoEdgeFace(neighbour_face, ioAffectedFaces))
  935. {
  936. // No, we need to recalculate its plane
  937. neighbour_face->CalculateNormalAndCentroid(mPositions.data());
  938. // Mark neighbour face as affected
  939. sMarkAffected(neighbour_face, ioAffectedFaces);
  940. }
  941. // Check if inFace now has only 2 edges left
  942. if (RemoveTwoEdgeFace(inFace, ioAffectedFaces))
  943. return; // Bail if face no longer exists
  944. // Restart loop
  945. recalculate_plane = true;
  946. removed = true;
  947. break;
  948. }
  949. // This edge is ok, go to the next edge
  950. edge = next_edge;
  951. neighbour_face = next_neighbour_face;
  952. } while (edge != inFace->mFirstEdge);
  953. } while (removed);
  954. // Recalculate plane?
  955. if (recalculate_plane)
  956. inFace->CalculateNormalAndCentroid(mPositions.data());
  957. }
  958. bool ConvexHullBuilder::RemoveTwoEdgeFace(Face *inFace, Faces &ioAffectedFaces) const
  959. {
  960. // Check if this face contains only 2 edges
  961. Edge *edge = inFace->mFirstEdge;
  962. Edge *next_edge = edge->mNextEdge;
  963. JPH_ASSERT(edge != next_edge); // 1 edge faces should not exist
  964. if (next_edge->mNextEdge == edge)
  965. {
  966. #ifdef JPH_CONVEX_BUILDER_DEBUG
  967. DrawWireFace(inFace, Color::sRed);
  968. DrawState();
  969. #endif
  970. // Schedule both neighbours for re-checking
  971. Edge *neighbour_edge = edge->mNeighbourEdge;
  972. Face *neighbour_face = neighbour_edge->mFace;
  973. Edge *next_neighbour_edge = next_edge->mNeighbourEdge;
  974. Face *next_neighbour_face = next_neighbour_edge->mFace;
  975. sMarkAffected(neighbour_face, ioAffectedFaces);
  976. sMarkAffected(next_neighbour_face, ioAffectedFaces);
  977. // Link my neighbours to each other
  978. neighbour_edge->mNeighbourEdge = next_neighbour_edge;
  979. next_neighbour_edge->mNeighbourEdge = neighbour_edge;
  980. // Unlink my edges
  981. edge->mNeighbourEdge = nullptr;
  982. next_edge->mNeighbourEdge = nullptr;
  983. // Mark this face as removed
  984. inFace->mRemoved = true;
  985. return true;
  986. }
  987. return false;
  988. }
  989. #ifdef JPH_ENABLE_ASSERTS
  990. void ConvexHullBuilder::DumpFace(const Face *inFace) const
  991. {
  992. Trace("f:0x%p", inFace);
  993. const Edge *e = inFace->mFirstEdge;
  994. do
  995. {
  996. Trace("\te:0x%p { i:%d e:0x%p f:0x%p }", e, e->mStartIdx, e->mNeighbourEdge, e->mNeighbourEdge->mFace);
  997. e = e->mNextEdge;
  998. } while (e != inFace->mFirstEdge);
  999. }
  1000. void ConvexHullBuilder::DumpFaces() const
  1001. {
  1002. Trace("Dump Faces:");
  1003. for (const Face *f : mFaces)
  1004. if (!f->mRemoved)
  1005. DumpFace(f);
  1006. }
  1007. void ConvexHullBuilder::ValidateFace(const Face *inFace) const
  1008. {
  1009. if (inFace->mRemoved)
  1010. {
  1011. const Edge *e = inFace->mFirstEdge;
  1012. if (e != nullptr)
  1013. do
  1014. {
  1015. JPH_ASSERT(e->mNeighbourEdge == nullptr);
  1016. e = e->mNextEdge;
  1017. } while (e != inFace->mFirstEdge);
  1018. }
  1019. else
  1020. {
  1021. int edge_count = 0;
  1022. const Edge *e = inFace->mFirstEdge;
  1023. do
  1024. {
  1025. // Count edge
  1026. ++edge_count;
  1027. // Validate that adjacent faces are all different
  1028. if (mFaces.size() > 2)
  1029. for (const Edge *other_edge = e->mNextEdge; other_edge != inFace->mFirstEdge; other_edge = other_edge->mNextEdge)
  1030. JPH_ASSERT(e->mNeighbourEdge->mFace != other_edge->mNeighbourEdge->mFace);
  1031. // Assert that the face is correct
  1032. JPH_ASSERT(e->mFace == inFace);
  1033. // Assert that we have a neighbour
  1034. const Edge *nb_edge = e->mNeighbourEdge;
  1035. JPH_ASSERT(nb_edge != nullptr);
  1036. if (nb_edge != nullptr)
  1037. {
  1038. // Assert that our neighbours edge points to us
  1039. JPH_ASSERT(nb_edge->mNeighbourEdge == e);
  1040. // Assert that it belongs to a different face
  1041. JPH_ASSERT(nb_edge->mFace != inFace);
  1042. // Assert that the next edge of the neighbour points to the same vertex as this edge's vertex
  1043. JPH_ASSERT(nb_edge->mNextEdge->mStartIdx == e->mStartIdx);
  1044. // Assert that my next edge points to the same vertex as my neighbours vertex
  1045. JPH_ASSERT(e->mNextEdge->mStartIdx == nb_edge->mStartIdx);
  1046. }
  1047. e = e->mNextEdge;
  1048. } while (e != inFace->mFirstEdge);
  1049. // Assert that we have 3 or more edges
  1050. JPH_ASSERT(edge_count >= 3);
  1051. }
  1052. }
  1053. void ConvexHullBuilder::ValidateFaces() const
  1054. {
  1055. for (const Face *f : mFaces)
  1056. ValidateFace(f);
  1057. }
  1058. #endif // JPH_ENABLE_ASSERTS
  1059. void ConvexHullBuilder::GetCenterOfMassAndVolume(Vec3 &outCenterOfMass, float &outVolume) const
  1060. {
  1061. // Fourth point is the average of all face centroids
  1062. Vec3 v4 = Vec3::sZero();
  1063. for (const Face *f : mFaces)
  1064. v4 += f->mCentroid;
  1065. v4 /= float(mFaces.size());
  1066. // Calculate mass and center of mass of this convex hull by summing all tetrahedrons
  1067. outVolume = 0.0f;
  1068. outCenterOfMass = Vec3::sZero();
  1069. for (const Face *f : mFaces)
  1070. {
  1071. // Get the first vertex that we'll use to create a triangle fan
  1072. Edge *e = f->mFirstEdge;
  1073. Vec3 v1 = mPositions[e->mStartIdx];
  1074. // Get the second vertex
  1075. e = e->mNextEdge;
  1076. Vec3 v2 = mPositions[e->mStartIdx];
  1077. for (e = e->mNextEdge; e != f->mFirstEdge; e = e->mNextEdge)
  1078. {
  1079. // Fetch the last point of the triangle
  1080. Vec3 v3 = mPositions[e->mStartIdx];
  1081. // Calculate center of mass and mass of this tetrahedron,
  1082. // see: https://en.wikipedia.org/wiki/Tetrahedron#Volume
  1083. 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
  1084. Vec3 center_of_mass_tetrahedron = v1 + v2 + v3 + v4; // Needs to be divided by 4, postpone this until the end of the loop
  1085. // Accumulate results
  1086. outVolume += volume_tetrahedron;
  1087. outCenterOfMass += volume_tetrahedron * center_of_mass_tetrahedron;
  1088. // Update v2 for next triangle
  1089. v2 = v3;
  1090. } while (e != f->mFirstEdge);
  1091. }
  1092. // Calculate center of mass, fall back to average point in case there is no volume (everything is on a plane in this case)
  1093. if (outVolume > FLT_EPSILON)
  1094. outCenterOfMass /= 4.0f * outVolume;
  1095. else
  1096. outCenterOfMass = v4;
  1097. outVolume /= 6.0f;
  1098. }
  1099. void ConvexHullBuilder::DetermineMaxError(Face *&outFaceWithMaxError, float &outMaxError, int &outMaxErrorPositionIdx, float &outCoplanarDistance) const
  1100. {
  1101. outCoplanarDistance = DetermineCoplanarDistance();
  1102. // This measures the distance from a polygon to the furthest point outside of the hull
  1103. float max_error = 0.0f;
  1104. Face *max_error_face = nullptr;
  1105. int max_error_point = -1;
  1106. for (int i = 0; i < (int)mPositions.size(); ++i)
  1107. {
  1108. Vec3 v = mPositions[i];
  1109. // This measures the closest edge from all faces to point v
  1110. // 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
  1111. // 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
  1112. // polygon.
  1113. float min_edge_dist_sq = FLT_MAX;
  1114. Face *min_edge_dist_face = nullptr;
  1115. for (Face *f : mFaces)
  1116. {
  1117. // Check if point is on or in front of plane
  1118. float normal_len = f->mNormal.Length();
  1119. JPH_ASSERT(normal_len > 0.0f);
  1120. float plane_dist = f->mNormal.Dot(v - f->mCentroid) / normal_len;
  1121. if (plane_dist > -outCoplanarDistance)
  1122. {
  1123. // Check distance to the edges of this face
  1124. float edge_dist_sq = GetDistanceToEdgeSq(v, f);
  1125. if (edge_dist_sq < min_edge_dist_sq)
  1126. {
  1127. min_edge_dist_sq = edge_dist_sq;
  1128. min_edge_dist_face = f;
  1129. }
  1130. // If the point is inside the polygon and the point is in front of the plane, measure the distance
  1131. if (edge_dist_sq == 0.0f && plane_dist > max_error)
  1132. {
  1133. max_error = plane_dist;
  1134. max_error_face = f;
  1135. max_error_point = i;
  1136. }
  1137. }
  1138. }
  1139. // If the minimum distance to an edge is further than our current max error, we use that as max error
  1140. float min_edge_dist = sqrt(min_edge_dist_sq);
  1141. if (min_edge_dist_face != nullptr && min_edge_dist > max_error)
  1142. {
  1143. max_error = min_edge_dist;
  1144. max_error_face = min_edge_dist_face;
  1145. max_error_point = i;
  1146. }
  1147. }
  1148. outFaceWithMaxError = max_error_face;
  1149. outMaxError = max_error;
  1150. outMaxErrorPositionIdx = max_error_point;
  1151. }
  1152. #ifdef JPH_CONVEX_BUILDER_DEBUG
  1153. void ConvexHullBuilder::DrawState(bool inDrawConflictList) const
  1154. {
  1155. // Draw origin
  1156. DebugRenderer::sInstance->DrawMarker(cDrawScale * mOffset, Color::sRed, 0.2f);
  1157. int face_idx = 0;
  1158. // Draw faces
  1159. for (const Face *f : mFaces)
  1160. if (!f->mRemoved)
  1161. {
  1162. Color iteration_color = Color::sGetDistinctColor(f->mIteration);
  1163. Color face_color = Color::sGetDistinctColor(face_idx++);
  1164. // First point
  1165. const Edge *e = f->mFirstEdge;
  1166. Vec3 p1 = cDrawScale * (mPositions[e->mStartIdx] + mOffset);
  1167. // Second point
  1168. e = e->mNextEdge;
  1169. Vec3 p2 = cDrawScale * (mPositions[e->mStartIdx] + mOffset);
  1170. // First line
  1171. DebugRenderer::sInstance->DrawLine(p1, p2, Color::sGrey);
  1172. do
  1173. {
  1174. // Third point
  1175. e = e->mNextEdge;
  1176. Vec3 p3 = cDrawScale * (mPositions[e->mStartIdx] + mOffset);
  1177. DebugRenderer::sInstance->DrawTriangle(p1, p2, p3, iteration_color);
  1178. DebugRenderer::sInstance->DrawLine(p2, p3, Color::sGrey);
  1179. p2 = p3;
  1180. }
  1181. while (e != f->mFirstEdge);
  1182. // Draw normal
  1183. Vec3 centroid = cDrawScale * (f->mCentroid + mOffset);
  1184. DebugRenderer::sInstance->DrawArrow(centroid, centroid + f->mNormal.NormalizedOr(Vec3::sZero()), face_color, 0.01f);
  1185. // Draw conflict list
  1186. if (inDrawConflictList)
  1187. for (int idx : f->mConflictList)
  1188. DebugRenderer::sInstance->DrawMarker(cDrawScale * (mOffset + mPositions[idx]), face_color, 0.05f);
  1189. }
  1190. // Offset to the right
  1191. mOffset += mDelta;
  1192. }
  1193. void ConvexHullBuilder::DrawWireFace(const Face *inFace, ColorArg inColor) const
  1194. {
  1195. const Edge *e = inFace->mFirstEdge;
  1196. Vec3 prev = cDrawScale * (mPositions[e->mStartIdx] + mOffset);
  1197. do
  1198. {
  1199. const Edge *next = e->mNextEdge;
  1200. Vec3 cur = cDrawScale * (mPositions[next->mStartIdx] + mOffset);
  1201. DebugRenderer::sInstance->DrawArrow(prev, cur, inColor, 0.01f);
  1202. DebugRenderer::sInstance->DrawText3D(prev, ConvertToString(e->mStartIdx), inColor);
  1203. e = next;
  1204. prev = cur;
  1205. } while (e != inFace->mFirstEdge);
  1206. }
  1207. void ConvexHullBuilder::DrawEdge(const Edge *inEdge, ColorArg inColor) const
  1208. {
  1209. Vec3 p1 = cDrawScale * (mPositions[inEdge->mStartIdx] + mOffset);
  1210. Vec3 p2 = cDrawScale * (mPositions[inEdge->mNextEdge->mStartIdx] + mOffset);
  1211. DebugRenderer::sInstance->DrawArrow(p1, p2, inColor, 0.01f);
  1212. }
  1213. #endif // JPH_CONVEX_BUILDER_DEBUG
  1214. #ifdef JPH_CONVEX_BUILDER_DUMP_SHAPE
  1215. void ConvexHullBuilder::DumpShape() const
  1216. {
  1217. static atomic<int> sShapeNo = 1;
  1218. int shape_no = sShapeNo++;
  1219. ofstream f;
  1220. f.open(StringFormat("dumped_shape%d.cpp", shape_no).c_str(), ofstream::out | ofstream::trunc);
  1221. if (!f.is_open())
  1222. return;
  1223. f << "{\n";
  1224. for (Vec3 v : mPositions)
  1225. f << StringFormat("\tVec3(%.9gf, %.9gf, %.9gf),\n", v.GetX(), v.GetY(), v.GetZ());
  1226. f << "},\n";
  1227. }
  1228. #endif // JPH_CONVEX_BUILDER_DUMP_SHAPE
  1229. JPH_NAMESPACE_END