Skip to content

Commit

Permalink
SimplifyQuadricDecimation fixes (#6163)
Browse files Browse the repository at this point in the history 
* SimplifyQuadricDecimation: Check for flipped triangles connected to both vertices of the removed edge

The previous code checks if any triangles connected to the removed vertex are flipped in the process of collapsing the edge. However, since the remaining vertex on the collapsed edge can be moved during the operation, triangles connected to this vertex can flip as well. This commit puts the flip-checking code into a loop so we can check both sets of triangles.

* SimplifyQuadricDecimation: Disallow creating very small triangles

This sets a threshold on how much any triangle can shrink in area during an edge collapse which should avoid creating degenerate triangles. The threshold is somewhat arbitrary, I found 0.001 worked well but it may make sense for this to become a function parameter.

* SimplifyQuadricDecimation: Prevent 'pinching' of triangles

When a vertex is connected to exactly 3 triangles and an edge of one of these triangles not connected to the vertex is collapsed, the remaining 2 triangles end up sharing the same vertices, with opposite normals. This creates a non-manifold location which is nearly always undesirable.
This commit adds a check to see if any of the triangles attached to the collapsing edge share an edge that is not connected to `vidx0` or `vidx1`

* SimplifyQuadricDecimation: Save some memory by clearing some vectors after they are no longer needed

* Update Changelog with SimplifyQuadricDecimation fixes

* SimplifyQuadricDecimation: Reduce memory requirements by removing need for vbars & costs maps

Reduces the memory footprint by 27% in my tests runs ~15% faster as we remove a lot of allocations.

* SimplifyQuadricDecimation: Simplify checking for pinched triangles

Changing from tracking edges (`unordered_set<Vector2i>`) to vertices (`unordered_map<int, int>) gives a good performance increase
  • Loading branch information
@bjudeworley
bjudeworley committed Aug 15, 2023
1 parent 75dccb4 commit a89031c
Show file tree
Hide file tree
Showing 2 changed files with 98 additions and 71 deletions.
1 change: 1 addition & 0 deletions CHANGELOG.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,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

Expand Down
168 changes: 97 additions & 71 deletions cpp/open3d/geometry/TriangleMeshSimplification.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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));
}
};
Expand All @@ -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();
Expand All @@ -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;
}

Expand Down Expand Up @@ -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] +
Expand Down

0 comments on commit a89031c

Please sign in to comment.