From b95a9f6e868ed9fd11a466c53e027fe000e3b8bc Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Mon, 12 Oct 2026 09:23:03 -0600 Subject: [PATCH] Fix shared boundary fan artifacts with multiple boundary loops (#2610) snap_boundary_to_mesh snapped whichever boundary loop libigl returned first, so a pre-existing loop (e.g. a cropped shaft) could be dragged onto the shared boundary. Snap only the seam vertices left by removing the shared surface, against all loops of the shared surface. extract_boundary_loop now picks the longest loop instead of an arbitrary one. Adds regression tests for both. --- Libs/Mesh/MeshUtils.cpp | 112 ++++++++++++++++++++---------- Testing/MeshTests/MeshTests.cpp | 116 ++++++++++++++++++++++++++++++++ 2 files changed, 193 insertions(+), 35 deletions(-) diff --git a/Libs/Mesh/MeshUtils.cpp b/Libs/Mesh/MeshUtils.cpp index 5ebe7ca8f4..4de45ae76c 100644 --- a/Libs/Mesh/MeshUtils.cpp +++ b/Libs/Mesh/MeshUtils.cpp @@ -313,7 +313,25 @@ Mesh MeshUtils::extract_boundary_loop(Mesh mesh) { throw std::runtime_error("Expected at least one boundary loop in the mesh"); } - auto loop = loops[0]; // copy so we can rotate it to a canonical start vertex + // A mesh may have more than one boundary loop (e.g. a hole in the shared surface). igl::boundary_loop + // returns them in arbitrary order, so use the longest one. + size_t longest = 0; + if (loops.size() > 1) { + double longest_length = -1.0; + for (size_t i = 0; i < loops.size(); i++) { + double length = 0.0; + for (size_t j = 0; j < loops[i].size(); j++) { + length += (V.row(loops[i][j]) - V.row(loops[i][(j + 1) % loops[i].size()])).norm(); + } + if (length > longest_length) { + longest_length = length; + longest = i; + } + } + SW_WARN("Mesh has {} boundary loops, extracting the longest", loops.size()); + } + + auto loop = loops[longest]; // copy so we can rotate it to a canonical start vertex // Rotate the loop so it always starts at a canonical vertex (the one with the // largest Y coordinate; lexicographic tiebreakers on Z then X). igl::boundary_loop @@ -377,9 +395,9 @@ Mesh MeshUtils::extract_boundary_loop(Mesh mesh) { */ //--------------------------------------------------------------------------- -static std::tuple rem_into_eigen_mesh(const std::vector& faces, - const Eigen::MatrixXd& src_V, - const Eigen::MatrixXi& src_F) { +// Also returns the list mapping each output vertex to its index in src_V +static std::tuple rem_into_eigen_mesh( + const std::vector& faces, const Eigen::MatrixXd& src_V, const Eigen::MatrixXi& src_F) { Eigen::MatrixXd V; Eigen::MatrixXi F; const std::unordered_set faces_set(faces.begin(), faces.end()); @@ -393,9 +411,10 @@ static std::tuple rem_into_eigen_mesh(const st } Eigen::VectorXi mapping; - igl::remove_unreferenced(src_V, F2, V, F, mapping); + Eigen::VectorXi new_to_old; + igl::remove_unreferenced(src_V, F2, V, F, mapping, new_to_old); - return std::make_tuple(V, F); + return std::make_tuple(V, F, new_to_old); } //--------------------------------------------------------------------------- @@ -425,18 +444,20 @@ static bool is_empty(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F) { retur //--------------------------------------------------------------------------- // Identifies faces from the source mesh that lie within tolerance distance of the other mesh surface. // Uses spatial queries to determine which faces are "shared" between the two meshes. -// Returns both the shared faces as one mesh and the remaining non-shared faces as another mesh. -static std::tuple find_shared_surface( - const Eigen::MatrixXd& src_V, const Eigen::MatrixXi& src_F, const Eigen::MatrixXd& other_V, - const Eigen::MatrixXi& other_F, double tol = 1e-3) { +// Returns both the shared faces as one mesh and the remaining non-shared faces as another mesh, along with the +// seam vertices: the vertices of the remaining mesh that were also used by a shared face. +static std::tuple> +find_shared_surface(const Eigen::MatrixXd& src_V, const Eigen::MatrixXi& src_F, const Eigen::MatrixXd& other_V, + const Eigen::MatrixXi& other_F, double tol = 1e-3) { Eigen::MatrixXd out_V; Eigen::MatrixXi out_F; Eigen::MatrixXd rem_V; Eigen::MatrixXi rem_F; + std::vector seam; // If either mesh is empty, there can be no shared surface if (is_empty(src_V, src_F) || is_empty(other_V, other_F)) { - return std::make_tuple(out_V, out_F, src_V, src_F); + return std::make_tuple(out_V, out_F, src_V, src_F, seam); } igl::AABB tree; @@ -470,36 +491,56 @@ static std::tuple in_shared(src_V.rows(), false); + for (const auto& f : shared_faces) { + for (int i = 0; i < 3; i++) { + in_shared[src_F(f, i)] = true; + } + } + for (int i = 0; i < rem_to_src.size(); i++) { + if (in_shared[rem_to_src(i)]) { + seam.push_back(i); + } + } + + return std::make_tuple(out_V, out_F, rem_V, rem_F, seam); } //--------------------------------------------------------------------------- -// Snaps boundary vertices of the source mesh to the closest points on the shared mesh boundary. +// Snaps the seam vertices of the source mesh to the closest points on the shared mesh boundary. // This ensures the boundaries are geometrically aligned, which is typically needed for mesh // stitching or boolean operations. -// Returns a copy of the source mesh with modified boundary vertex positions. +// Only the seam vertices (those left on the boundary by removing the shared surface) are moved, so any other +// boundary loops of the source mesh (e.g. a cropped bone shaft) are left untouched. All boundary loops of the +// shared mesh are considered, since the shared surface may have holes or consist of several patches. +// Returns a copy of the source mesh with modified seam vertex positions. static void snap_boundary_to_mesh(const Eigen::MatrixXd& src_V, const Eigen::MatrixXi& src_F, - const Eigen::MatrixXd& shared_V, const Eigen::MatrixXi& shared_F, - Eigen::MatrixXd& out_V, Eigen::MatrixXi& out_F) { - std::vector > src_loops, shared_loops; - igl::boundary_loop(src_F, src_loops); + const std::vector& seam, const Eigen::MatrixXd& shared_V, + const Eigen::MatrixXi& shared_F, Eigen::MatrixXd& out_V, Eigen::MatrixXi& out_F) { + std::vector > shared_loops; igl::boundary_loop(shared_F, shared_loops); - if (src_loops.size() < 1 || shared_loops.size() < 1) { - SW_DEBUG("src_loops size: {}, shared_loops size: {}", src_loops.size(), shared_loops.size()); - throw std::runtime_error("Expected at least one boundary loop in each mesh"); + if (shared_loops.size() < 1) { + throw std::runtime_error("Expected at least one boundary loop in the shared surface"); } - const auto& src_loop = src_loops[0]; - const auto& shared_loop = shared_loops[0]; + size_t num_edges = 0; + for (const auto& shared_loop : shared_loops) { + num_edges += shared_loop.size(); + } Eigen::MatrixXi shared_F_boundary; - shared_F_boundary.resize(shared_loop.size(), 3); - for (int i = 0; i < shared_loop.size(); i++) { - const int v0 = shared_loop[i]; - const int v1 = shared_loop[(i + 1) % shared_loop.size()]; - shared_F_boundary.row(i) = Eigen::Vector3i{v0, v1, v1}; + shared_F_boundary.resize(num_edges, 3); + int next_idx = 0; + for (const auto& shared_loop : shared_loops) { + for (int i = 0; i < shared_loop.size(); i++) { + const int v0 = shared_loop[i]; + const int v1 = shared_loop[(i + 1) % shared_loop.size()]; + shared_F_boundary.row(next_idx++) = Eigen::Vector3i{v0, v1, v1}; + } } out_V = src_V; @@ -508,15 +549,15 @@ static void snap_boundary_to_mesh(const Eigen::MatrixXd& src_V, const Eigen::Mat igl::AABB tree; tree.init(shared_V, shared_F_boundary); - for (int i = 0; i < src_loop.size(); i++) { + for (const auto& v : seam) { Eigen::VectorXd sqrD; Eigen::VectorXi I; Eigen::MatrixXd C; Eigen::MatrixXd P(1, 3); - P.row(0) = src_V.row(src_loop[i]); + P.row(0) = src_V.row(v); tree.squared_distance(shared_V, shared_F_boundary, P, sqrD, I, C); - out_V.row(src_loop[i]) = C.row(0); + out_V.row(v) = C.row(0); } } @@ -538,8 +579,9 @@ std::array MeshUtils::shared_boundary_extractor(const Mesh& mesh_l, con Eigen::MatrixXd shared_V_l, shared_V_r, rem_V_l, rem_V_r; Eigen::MatrixXi shared_F_l, shared_F_r, rem_F_l, rem_F_r; - std::tie(shared_V_l, shared_F_l, rem_V_l, rem_F_l) = find_shared_surface(V_l, F_l, V_r, F_r, tol); - std::tie(shared_V_r, shared_F_r, rem_V_r, rem_F_r) = find_shared_surface(V_r, F_r, V_l, F_l, tol); + std::vector seam_l; + std::tie(shared_V_l, shared_F_l, rem_V_l, rem_F_l, seam_l) = find_shared_surface(V_l, F_l, V_r, F_r, tol); + std::tie(shared_V_r, shared_F_r, rem_V_r, rem_F_r, std::ignore) = find_shared_surface(V_r, F_r, V_l, F_l, tol); if (is_empty(shared_V_l, shared_F_l) || is_empty(shared_V_r, shared_F_r) || is_empty(rem_V_l, rem_F_l) || is_empty(rem_V_r, rem_F_r)) { @@ -551,7 +593,7 @@ std::array MeshUtils::shared_boundary_extractor(const Mesh& mesh_l, con Eigen::MatrixXd bridge_V; Eigen::MatrixXi bridge_F; - snap_boundary_to_mesh(rem_V_l, rem_F_l, shared_V_r, shared_F_r, bridge_V, bridge_F); + snap_boundary_to_mesh(rem_V_l, rem_F_l, seam_l, shared_V_r, shared_F_r, bridge_V, bridge_F); Mesh out_l(bridge_V, bridge_F); Mesh out_r(rem_V_r, rem_F_r); diff --git a/Testing/MeshTests/MeshTests.cpp b/Testing/MeshTests/MeshTests.cpp index 59ee0dd004..2ab85547e5 100644 --- a/Testing/MeshTests/MeshTests.cpp +++ b/Testing/MeshTests/MeshTests.cpp @@ -1,3 +1,4 @@ +#include #include #include #include @@ -1195,6 +1196,121 @@ TEST(MeshTests, boundaryLoopExtractor) { ASSERT_TRUE(loop == ground_truth); } +// Returns a copy of the mesh with the faces around the given vertex removed, leaving a hole +static Mesh punch_hole(const Mesh& mesh, int center) { + Eigen::MatrixXd V = mesh.points(); + Eigen::MatrixXi F = mesh.faces(); + Eigen::MatrixXi F2(F.rows(), 3); + int count = 0; + for (int f = 0; f < F.rows(); f++) { + if (F(f, 0) != center && F(f, 1) != center && F(f, 2) != center) { + F2.row(count++) = F.row(f); + } + } + F2.conservativeResize(count, 3); + return Mesh(V, F2); +} + +static double max_edge_length(const Mesh& mesh) { + Eigen::MatrixXd V = mesh.points(); + Eigen::MatrixXi F = mesh.faces(); + double max_length = 0.0; + for (int f = 0; f < F.rows(); f++) { + for (int i = 0; i < 3; i++) { + max_length = std::max(max_length, (V.row(F(f, i)) - V.row(F(f, (i + 1) % 3))).norm()); + } + } + return max_length; +} + +// A mesh that already has a boundary loop (e.g. a bone cropped at the shaft) must not have that loop dragged +// onto the shared boundary (issue #2610) +TEST(MeshTests, sharedBoundaryExtractorMultipleLoops) { + Mesh left(std::string(TEST_DATA_DIR) + "/shared_boundary/00_l.vtk"); + Mesh right(std::string(TEST_DATA_DIR) + "/shared_boundary/00_r.vtk"); + double tol = 0.001; + + Eigen::MatrixXd V = left.points(); + Eigen::MatrixXi F = left.faces(); + Eigen::VectorXd sqrD; + Eigen::VectorXi I; + Eigen::MatrixXd C; + igl::point_mesh_squared_distance(V, right.points(), right.faces(), sqrD, I, C); + + // The lowest indexed vertex anywhere near the shared surface. A hole whose rim has a lower indexed vertex + // is returned by igl::boundary_loop ahead of the loop left by removing the shared surface. + int first_near = V.rows(); + for (int i = 0; i < V.rows(); i++) { + if (sqrD(i) <= tol) { + first_near = i; + break; + } + } + + // Find a face far from the shared surface with a vertex below first_near, and punch a hole next to it + int center = -1; + for (int f = 0; f < F.rows() && center < 0; f++) { + for (int i = 0; i < 3; i++) { + int rim = F(f, i); + int candidate = F(f, (i + 1) % 3); + if (rim < first_near && sqrD(rim) > 0.25 * sqrD.maxCoeff() && sqrD(candidate) > 0.25 * sqrD.maxCoeff()) { + center = candidate; + break; + } + } + } + ASSERT_GE(center, 0); + + Mesh left_open = punch_hole(left, center); + std::array output = MeshUtils::shared_boundary_extractor(left_open, right, tol); + + // snapping only moves vertices a short distance, so no edge should end up spanning the mesh + ASSERT_LT(max_edge_length(output[0]), 3.0 * max_edge_length(left_open)); +} + +// With more than one boundary loop, the longest should be extracted rather than an arbitrary one +TEST(MeshTests, boundaryLoopExtractorMultipleLoops) { + Mesh ground_truth(std::string(TEST_DATA_DIR) + "/shared_boundary/00_out_c.vtp"); + Mesh mesh(std::string(TEST_DATA_DIR) + "/shared_boundary/00_out_s.vtk"); + + // punch a hole in the interior of the mesh + Eigen::MatrixXd V; + Eigen::MatrixXi F; + std::vector> loops; + for (int center = 0; center < mesh.numPoints(); center++) { + Mesh holed = punch_hole(mesh, center); + V = holed.points(); + F = holed.faces(); + loops.clear(); + igl::boundary_loop(F, loops); + if (loops.size() == 2) { + break; + } + } + ASSERT_EQ(loops.size(), 2); + + // igl::boundary_loop returns the loop with the lowest indexed vertex first, so swap vertex 0 with a vertex on + // the rim of the hole to make sure the hole is returned ahead of the outer loop + const int rim = loops[0].size() < loops[1].size() ? loops[0][0] : loops[1][0]; + V.row(0).swap(V.row(rim)); + for (int f = 0; f < F.rows(); f++) { + for (int i = 0; i < 3; i++) { + if (F(f, i) == 0) { + F(f, i) = rim; + } else if (F(f, i) == rim) { + F(f, i) = 0; + } + } + } + loops.clear(); + igl::boundary_loop(F, loops); + ASSERT_EQ(loops.size(), 2); + ASSERT_LT(loops[0].size(), loops[1].size()); + + Mesh loop = MeshUtils::extract_boundary_loop(Mesh(V, F)); + ASSERT_TRUE(loop == ground_truth); +} + TEST(MeshTests, thicknessTest) { Mesh sphere(std::string(TEST_DATA_DIR) + "/thickness/sphere.vtk"); Image ct(std::string(TEST_DATA_DIR) + "/thickness/ct.nrrd");