SimplifyQuadricDecimation fixes

CHANGELOG.md

@@ -12,6 +12,7 @@
 * Fix raycasting scene: Allow setting of number of threads that are used for building a raycasting scene
 * Fix Python bindings for CUDA device synchronization, voxel grid saving (PR #5425)
 * Support msgpack versions without cmake
+* Fix some bad triangle generation in TriangleMesh::SimplifyQuadricDecimation
 
 ## 0.13
 

cpp/open3d/geometry/TriangleMeshSimplification.cpp

@@ -311,51 +311,55 @@ std::shared_ptr<TriangleMesh> TriangleMesh::SimplifyQuadricDecimation(
  AddPerpPlaneQuadric(tria(2), tria(0), tria(1), area);
  }
 
+ // Clear the unused vectors to save some memory
+ triangle_areas.clear();
+ triangle_planes.clear();
+ edge_triangle_count.clear();
+
  // Get valid edges and compute cost
- // Note: We could also select all vertex pairs as edges with dist < eps
- std::unordered_map<Eigen::Vector2i, Eigen::Vector3d,
- utility::hash_eigen<Eigen::Vector2i>>
- vbars;
- std::unordered_map<Eigen::Vector2i, double,
- utility::hash_eigen<Eigen::Vector2i>>
- costs;
  auto CostEdgeComp = [](const CostEdge& a, const CostEdge& b) {
  return std::get<0>(a) > std::get<0>(b);
  };
  std::priority_queue<CostEdge, std::vector<CostEdge>, decltype(CostEdgeComp)>
  queue(CostEdgeComp);
 
- auto AddEdge = [&](int vidx0, int vidx1, bool update) {
- int min = std::min(vidx0, vidx1);
- int max = std::max(vidx0, vidx1);
- Eigen::Vector2i edge(min, max);
- if (update || vbars.count(edge) == 0) {
- const Quadric& Q0 = Qs[min];
- const Quadric& Q1 = Qs[max];
- Quadric Qbar = Q0 + Q1;
- double cost;
- Eigen::Vector3d vbar;
- if (Qbar.IsInvertible()) {
- vbar = Qbar.Minimum();
- cost = Qbar.Eval(vbar);
+ auto compute_cost_vbar = [&](Eigen::Vector2i e) {
+ const Quadric& Q0 = Qs[e(0)];
+ const Quadric& Q1 = Qs[e(1)];
+ const Quadric Qbar = Q0 + Q1;
+ double cost;
+ Eigen::Vector3d vbar;
+ if (Qbar.IsInvertible()) {
+ vbar = Qbar.Minimum();
+ cost = Qbar.Eval(vbar);
+ } else {
+ const Eigen::Vector3d& v0 = mesh->vertices_[e(0)];
+ const Eigen::Vector3d& v1 = mesh->vertices_[e(1)];
+ const Eigen::Vector3d vmid = (v0 + v1) / 2;
+ const double cost0 = Qbar.Eval(v0);
+ const double cost1 = Qbar.Eval(v1);
+ const double costmid = Qbar.Eval(vmid);
+ cost = std::min(cost0, std::min(cost1, costmid));
+ if (cost == costmid) {
+ vbar = vmid;
+ } else if (cost == cost0) {
+ vbar = v0;
  } else {
- const Eigen::Vector3d& v0 = mesh->vertices_[vidx0];
- const Eigen::Vector3d& v1 = mesh->vertices_[vidx1];
- Eigen::Vector3d vmid = (v0 + v1) / 2;
- double cost0 = Qbar.Eval(v0);
- double cost1 = Qbar.Eval(v1);
- double costmid = Qbar.Eval(vmid);
- cost = std::min(cost0, std::min(cost1, costmid));
- if (cost == costmid) {
- vbar = vmid;
- } else if (cost == cost0) {
- vbar = v0;
- } else {
- vbar = v1;
- }
+ vbar = v1;
  }
- vbars[edge] = vbar;
- costs[edge] = cost;
+ }
+ return std::make_pair(cost, vbar);
+ };
+
+ std::unordered_set<Eigen::Vector2i, utility::hash_eigen<Eigen::Vector2i>>
+ added_edges;
+ auto AddEdge = [&](int vidx0, int vidx1, bool update) {
+ const int min = std::min(vidx0, vidx1);
+ const int max = std::max(vidx0, vidx1);
+ const Eigen::Vector2i edge(min, max);
+ if (update || added_edges.count(edge) == 0) {
+ const auto cost = compute_cost_vbar(edge).first;
+ added_edges.insert(edge);
  queue.push(CostEdge(cost, min, max));
  }
  };
@@ -366,6 +370,7 @@ std::shared_ptr<TriangleMesh> TriangleMesh::SimplifyQuadricDecimation(
  AddEdge(triangle(1), triangle(2), false);
  AddEdge(triangle(2), triangle(0), false);
  }
+ added_edges.clear();
 
  // perform incremental edge collapse
  bool has_vert_normal = HasVertexNormals();
@@ -384,50 +389,71 @@ std::shared_ptr<TriangleMesh> TriangleMesh::SimplifyQuadricDecimation(
 
  // test if the edge has been updated (reinserted into queue)
  Eigen::Vector2i edge(vidx0, vidx1);
+ const auto cost_vbar = compute_cost_vbar(edge);
+ const Eigen::Vector3d vbar = cost_vbar.second;
  bool valid = !vertices_deleted[vidx0] && !vertices_deleted[vidx1] &&
- cost == costs[edge];
+ cost == cost_vbar.first;
  if (!valid) {
  continue;
  }
 
- // avoid flip of triangle normal
- bool flipped = false;
- for (int tidx : vert_to_triangles[vidx1]) {
- if (triangles_deleted[tidx]) {
- continue;
- }
-
- const Eigen::Vector3i& tria = mesh->triangles_[tidx];
- bool has_vidx0 =
- vidx0 == tria(0) || vidx0 == tria(1) || vidx0 == tria(2);
- bool has_vidx1 =
- vidx1 == tria(0) || vidx1 == tria(1) || vidx1 == tria(2);
- if (has_vidx0 && has_vidx1) {
- continue;
- }
+ // avoid flip of triangle normal and creation of degenerate triangles
+ bool creates_invalid_triangle = false;
+ const double degenerate_ratio_threshold = 0.001;
+ std::unordered_map<int, int> edges{};
+ for (int vidx : {vidx1, vidx0}) {
+ for (int tidx : vert_to_triangles[vidx]) {
+ if (triangles_deleted[tidx]) {
+ continue;
+ }
 
- Eigen::Vector3d vert0 = mesh->vertices_[tria(0)];
- Eigen::Vector3d vert1 = mesh->vertices_[tria(1)];
- Eigen::Vector3d vert2 = mesh->vertices_[tria(2)];
- Eigen::Vector3d norm_before = (vert1 - vert0).cross(vert2 - vert0);
- norm_before /= norm_before.norm();
+ const Eigen::Vector3i& tria = mesh->triangles_[tidx];
+ const bool has_vidx0 = vidx0 == tria(0) || vidx0 == tria(1) ||
+ vidx0 == tria(2);
+ const bool has_vidx1 = vidx1 == tria(0) || vidx1 == tria(1) ||
+ vidx1 == tria(2);
+ if (has_vidx0 && has_vidx1) {
+ continue;
+ }
 
- if (vidx1 == tria(0)) {
- vert0 = vbars[edge];
- } else if (vidx1 == tria(1)) {
- vert1 = vbars[edge];
- } else if (vidx1 == tria(2)) {
- vert2 = vbars[edge];
- }
+ Eigen::Vector3d verts[3] = {mesh->vertices_[tria(0)],
+ mesh->vertices_[tria(1)],
+ mesh->vertices_[tria(2)]};
+ Eigen::Vector3d norm_before =
+ (verts[1] - verts[0]).cross(verts[2] - verts[0]);
+ const double area_before = 0.5 * norm_before.norm();
+ norm_before /= norm_before.norm();
+
+ for (auto i = 0; i < 3; ++i) {
+ if (tria(i) == vidx) {
+ verts[i] = vbar;
+ continue;
+ }
+ auto& vert_count = edges[tria(i)];
+ creates_invalid_triangle |= vert_count >= 2;
+ vert_count += 1;
+ }
 
- Eigen::Vector3d norm_after = (vert1 - vert0).cross(vert2 - vert0);
- norm_after /= norm_after.norm();
- if (norm_before.dot(norm_after) < 0) {
- flipped = true;
- break;
+ Eigen::Vector3d norm_after =
+ (verts[1] - verts[0]).cross(verts[2] - verts[0]);
+ const double area_after = 0.5 * norm_after.norm();
+ norm_after /= norm_after.norm();
+ // Disallow flipping the triangle normal
+ creates_invalid_triangle |= norm_before.dot(norm_after) < 0;
+ // Disallow creating very small triangles (possibly degenerate)
+ creates_invalid_triangle |=
+ area_after < degenerate_ratio_threshold * area_before;
+
+ if (creates_invalid_triangle) {
+ // Goto is the only way to jump out of two loops without
+ // multiple redundant if()'s. Yes, it can lead to spagetti
+ // code if abused but we're doing a very short jump here
+ goto end_flip_loop;
+ }
  }
  }
- if (flipped) {
+ end_flip_loop:
+ if (creates_invalid_triangle) {
  continue;
  }
 
@@ -460,7 +486,7 @@ std::shared_ptr<TriangleMesh> TriangleMesh::SimplifyQuadricDecimation(
  }
 
  // update vertex vidx0 to vbar
- mesh->vertices_[vidx0] = vbars[edge];
+ mesh->vertices_[vidx0] = vbar;
  Qs[vidx0] += Qs[vidx1];
  if (has_vert_normal) {
  mesh->vertex_normals_[vidx0] = 0.5 * (mesh->vertex_normals_[vidx0] +