From 89d0934f71ab8593867d037ca3bfe96cd3491fa9 Mon Sep 17 00:00:00 2001 From: Per Melin Date: Thu, 11 Apr 2024 16:03:38 +0200 Subject: [PATCH] Delaunay3D: Improve triangulation --- core/math/delaunay_3d.h | 77 ++++++++++++++++++++--------------------- 1 file changed, 38 insertions(+), 39 deletions(-) diff --git a/core/math/delaunay_3d.h b/core/math/delaunay_3d.h index 45571fb57085..25bd4e8d892c 100644 --- a/core/math/delaunay_3d.h +++ b/core/math/delaunay_3d.h @@ -46,7 +46,8 @@ class Delaunay3D { struct Simplex; enum { - ACCEL_GRID_SIZE = 16 + ACCEL_GRID_SIZE = 16, + QUANTIZATION_MAX = 1 << 16 // A power of two smaller than the 23 bit significand of a float. }; struct GridPos { Vector3i pos; @@ -173,38 +174,25 @@ class Delaunay3D { R128 radius2 = rel2_x * rel2_x + rel2_y * rel2_y + rel2_z * rel2_z; - return radius2 < (p_simplex.circum_r2 - R128(0.00001)); + return radius2 < (p_simplex.circum_r2 - R128(0.0000000001)); + // When this tolerance is too big, it can result in overlapping simplices. + // When it's too small, large amounts of planar simplices are created. } static bool simplex_is_coplanar(const Vector3 *p_points, const Simplex &p_simplex) { - Plane p(p_points[p_simplex.points[0]], p_points[p_simplex.points[1]], p_points[p_simplex.points[2]]); - if (ABS(p.distance_to(p_points[p_simplex.points[3]])) < CMP_EPSILON) { - return true; + // Checking every possible distance like this is overkill, but only checking + // one is not enough. If the simplex is almost planar then the vectors p1-p2 + // and p1-p3 can be practically collinear, which makes Plane unreliable. + for (uint32_t i = 0; i < 4; i++) { + Plane p(p_points[p_simplex.points[i]], p_points[p_simplex.points[(i + 1) % 4]], p_points[p_simplex.points[(i + 2) % 4]]); + // This tolerance should not be smaller than the one used with + // Plane::has_point() when creating the LightmapGI probe BSP tree. + if (ABS(p.distance_to(p_points[p_simplex.points[(i + 3) % 4]])) < 0.001) { + return true; + } } - Projection cm; - - cm.columns[0][0] = p_points[p_simplex.points[0]].x; - cm.columns[0][1] = p_points[p_simplex.points[1]].x; - cm.columns[0][2] = p_points[p_simplex.points[2]].x; - cm.columns[0][3] = p_points[p_simplex.points[3]].x; - - cm.columns[1][0] = p_points[p_simplex.points[0]].y; - cm.columns[1][1] = p_points[p_simplex.points[1]].y; - cm.columns[1][2] = p_points[p_simplex.points[2]].y; - cm.columns[1][3] = p_points[p_simplex.points[3]].y; - - cm.columns[2][0] = p_points[p_simplex.points[0]].z; - cm.columns[2][1] = p_points[p_simplex.points[1]].z; - cm.columns[2][2] = p_points[p_simplex.points[2]].z; - cm.columns[2][3] = p_points[p_simplex.points[3]].z; - - cm.columns[3][0] = 1.0; - cm.columns[3][1] = 1.0; - cm.columns[3][2] = 1.0; - cm.columns[3][3] = 1.0; - - return ABS(cm.determinant()) <= CMP_EPSILON; + return false; } public: @@ -215,9 +203,10 @@ class Delaunay3D { static Vector tetrahedralize(const Vector &p_points) { uint32_t point_count = p_points.size(); Vector3 *points = (Vector3 *)memalloc(sizeof(Vector3) * (point_count + 4)); + const Vector3 *src_points = p_points.ptr(); + Vector3 proportions; { - const Vector3 *src_points = p_points.ptr(); AABB rect; for (uint32_t i = 0; i < point_count; i++) { Vector3 point = src_points[i]; @@ -226,17 +215,25 @@ class Delaunay3D { } else { rect.expand_to(point); } - points[i] = point; } + real_t longest_axis = rect.size[rect.get_longest_axis_index()]; + proportions = Vector3(longest_axis, longest_axis, longest_axis) / rect.size; + for (uint32_t i = 0; i < point_count; i++) { - points[i] = (points[i] - rect.position) / rect.size; + // Scale points to the unit cube to better utilize R128 precision + // and quantize to stabilize triangulation over a wide range of + // distances. + points[i] = Vector3(Vector3i((src_points[i] - rect.position) / longest_axis * QUANTIZATION_MAX)) / QUANTIZATION_MAX; } - const real_t delta_max = Math::sqrt(2.0) * 20.0; + const real_t delta_max = Math::sqrt(2.0) * 100.0; Vector3 center = Vector3(0.5, 0.5, 0.5); - // any simplex that contains everything is good + // The larger the root simplex is, the more likely it is that the + // triangulation is convex. If it's not absolutely huge, there can + // be missing simplices that are not created for the outermost faces + // of the point cloud if the point density is very low there. points[point_count + 0] = center + Vector3(0, 1, 0) * delta_max; points[point_count + 1] = center + Vector3(0, -1, 1) * delta_max; points[point_count + 2] = center + Vector3(1, -1, -1) * delta_max; @@ -271,7 +268,7 @@ class Delaunay3D { for (uint32_t i = 0; i < point_count; i++) { bool unique = true; for (uint32_t j = i + 1; j < point_count; j++) { - if (points[i].is_equal_approx(points[j])) { + if (points[i] == points[j]) { unique = false; break; } @@ -280,7 +277,7 @@ class Delaunay3D { continue; } - Vector3i grid_pos = Vector3i(points[i] * ACCEL_GRID_SIZE); + Vector3i grid_pos = Vector3i(points[i] * proportions * ACCEL_GRID_SIZE); grid_pos = grid_pos.clamp(Vector3i(), Vector3i(ACCEL_GRID_SIZE - 1, ACCEL_GRID_SIZE - 1, ACCEL_GRID_SIZE - 1)); for (List::Element *E = acceleration_grid[grid_pos.x][grid_pos.y][grid_pos.z].front(); E;) { @@ -300,6 +297,9 @@ class Delaunay3D { Triangle t = Triangle(simplex->points[triangle_order[k][0]], simplex->points[triangle_order[k][1]], simplex->points[triangle_order[k][2]]); uint32_t *p = triangles_inserted.lookup_ptr(t); if (p) { + // This Delaunay implementation uses the Bowyer-Watson algorithm. + // The rule is that you don't reuse any triangles that were + // shared by any of the retriangulated simplices. triangles[*p].bad = true; } else { triangles_inserted.insert(t, triangles.size()); @@ -307,7 +307,6 @@ class Delaunay3D { } } - //remove simplex and continue simplex_list.erase(simplex->SE); for (const GridPos &gp : simplex->grid_positions) { @@ -334,8 +333,8 @@ class Delaunay3D { const real_t radius2 = Math::sqrt(double(new_simplex->circum_r2)) + 0.0001; Vector3 extents = Vector3(radius2, radius2, radius2); - Vector3i from = Vector3i((center - extents) * ACCEL_GRID_SIZE); - Vector3i to = Vector3i((center + extents) * ACCEL_GRID_SIZE); + Vector3i from = Vector3i((center - extents) * proportions * ACCEL_GRID_SIZE); + Vector3i to = Vector3i((center + extents) * proportions * ACCEL_GRID_SIZE); from = from.clamp(Vector3i(), Vector3i(ACCEL_GRID_SIZE - 1, ACCEL_GRID_SIZE - 1, ACCEL_GRID_SIZE - 1)); to = to.clamp(Vector3i(), Vector3i(ACCEL_GRID_SIZE - 1, ACCEL_GRID_SIZE - 1, ACCEL_GRID_SIZE - 1)); @@ -370,7 +369,7 @@ class Delaunay3D { break; } } - if (invalid || simplex_is_coplanar(points, *simplex)) { + if (invalid || simplex_is_coplanar(src_points, *simplex)) { memdelete(simplex); continue; }