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