From e94be55f46d07febfca46f892a5b01cfa9093cdd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Mon, 22 Jan 2024 18:04:52 +0100 Subject: [PATCH] add prototype for a function to snap vertices of 2 meshes onto each other faces --- .../Polygon_mesh_processing/CMakeLists.txt | 1 + .../internal/surface_snapping.h | 268 ++++++++++++++++++ .../Polygon_mesh_processing/CMakeLists.txt | 1 + .../snap_and_corefinement_union.cpp | 44 +++ 4 files changed, 314 insertions(+) create mode 100644 Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/surface_snapping.h create mode 100644 Polygon_mesh_processing/test/Polygon_mesh_processing/snap_and_corefinement_union.cpp diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt index 455e1ce2be7a..8c32d7cad8ff 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt @@ -34,6 +34,7 @@ create_single_source_cgal_program("corefinement_SM.cpp") create_single_source_cgal_program("corefinement_consecutive_bool_op.cpp") create_single_source_cgal_program("corefinement_difference_remeshed.cpp") create_single_source_cgal_program("corefinement_mesh_union.cpp") +create_single_source_cgal_program("snap_and_corefinement_union.cpp") create_single_source_cgal_program("corefinement_mesh_union_progress.cpp") create_single_source_cgal_program("corefinement_mesh_union_and_intersection.cpp") create_single_source_cgal_program("corefinement_mesh_union_with_attributes.cpp") diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/surface_snapping.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/surface_snapping.h new file mode 100644 index 000000000000..2e1d3beebace --- /dev/null +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/surface_snapping.h @@ -0,0 +1,268 @@ +// Copyright (c) 2016 GeometryFactory (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// +// Author(s) : Sebastien Loriot + +#ifndef CGAL_POLYGON_MESH_PROCESSING_INTERNAL_SURFACE_SNAPPING_H +#define CGAL_POLYGON_MESH_PROCESSING_INTERNAL_SURFACE_SNAPPING_H + +#include + +#include +#include +#include + +#include + + +namespace CGAL { +namespace Polygon_mesh_processing { +namespace experimental { + +template +void surface_snapping_impl_tm1_on_tm2(const AABB_tree_& tm1_tree, + TriangleMesh& tm1, TriangleMesh& tm2, + VPM vpm1, VPM vpm2, + MatchingMap matching_vertices_1, MatchingMap matching_vertices_2, + double eps2) +{ + using BGT = boost::graph_traits; + using vertex_descriptor = typename BGT::vertex_descriptor; + using face_descriptor = typename BGT::face_descriptor; + using halfedge_descriptor = typename BGT::halfedge_descriptor; + + auto split_edge_and_incident_faces = + [](halfedge_descriptor h, const typename K::Point_3& pt, TriangleMesh& tm, VPM vpm) + { + halfedge_descriptor res = Euler::split_edge(h, tm); + put(vpm, target(res, tm), pt); + + if(!is_border(res, tm)) + Euler::split_face(res, next(h, tm), tm); + + halfedge_descriptor opp_h = opposite(h, tm); + if(!is_border(opp_h, tm)) + Euler::split_face(opp_h, next(next(opp_h, tm), tm), tm); + + return res; + }; + + // TODO PARALLEL FOR + for (vertex_descriptor vd : vertices(tm2)) + { + if (get(matching_vertices_2, vd)) continue; // already on tm1 + typename K::Sphere_3 epsilon_sphere(get(vpm2, vd), eps2); + std::vector closed_faces; + tm1_tree.all_intersected_primitives(epsilon_sphere, std::back_inserter(closed_faces)); + + switch(closed_faces.size()) + { + case 0: + break; + case 1: + { + // check if close to border edges + face_descriptor f1 = closed_faces[0]; + halfedge_descriptor h1 = halfedge(f1, tm1); + std::array f1_hedges; + std::vector closed_hedges; + for (int i=0; i<3; ++i) + { + f1_hedges[i]=h1; + if (face(opposite(h1, tm1), tm1) == BGT::null_face()) + { + if (do_intersect(typename K::Segment_3(get(vpm1,source(h1,tm1)), + get(vpm1,target(h1,tm1))), + epsilon_sphere)) + { + closed_hedges.push_back(i); + } + } + h1 = next(h1, tm1); + } + switch(closed_hedges.size()) + { + case 0: + { + halfedge_descriptor hnew = Euler::add_center_vertex(h1,tm1); + put(vpm1, target(hnew, tm1), get(vpm2, vd)); + put(matching_vertices_1, target(hnew, tm1), true); + put(matching_vertices_2, vd, true); + break; + } + case 1: + { + h1 = f1_hedges[closed_hedges[0]]; + // try to snap on h1 edge points + if (compare_squared_distance(get(vpm1, target(h1, tm1)), get(vpm2, vd), eps2)!=LARGER) + { + put(vpm2, vd, get(vpm1, target(h1, tm1))); + put(matching_vertices_1, target(h1, tm1), true); + put(matching_vertices_2, vd, true); + break; + } + if (compare_squared_distance(get(vpm1, source(h1, tm1)), get(vpm2, vd), eps2)!=LARGER) + { + put(vpm2, vd, get(vpm1, source(h1, tm1))); + put(matching_vertices_1, source(h1, tm1), true); + put(matching_vertices_2, vd, true); + break; + } + // split and snap onto the edge + halfedge_descriptor hnew = + split_edge_and_incident_faces(h1, get(vpm2, vd), tm1, vpm1); + put(matching_vertices_1, target(hnew, tm1), true); + put(matching_vertices_2, vd, true); + break; + } + case 2: + { + // snap onto the common vertex + h1=next(f1_hedges[closed_hedges[0]+closed_hedges[1]-3], tm1); + put(vpm2, vd, get(vpm1, target(h1, tm1))); + put(matching_vertices_1, target(h1, tm1), true); + put(matching_vertices_2, vd, true); + break; + } + default: + //snap to any of the face vertex + put(vpm2, vd, get(vpm1,target(h1, tm1))); + put(matching_vertices_1, target(h1, tm1), true); + put(matching_vertices_2, vd, true); + } + break; + } + case 2: + { + // TODO: in case some edges are boundary edges, we might snap to a vertex + // look for a common edge + face_descriptor f1 = closed_faces[0], f2 = closed_faces[1]; + halfedge_descriptor h1 = halfedge(f1, tm1); + halfedge_descriptor h_to_split = BGT::null_halfedge(); + for (int i=0; i<3; ++i) + { + if (face(opposite(h1, tm1), tm1) == f2) + { + h_to_split=h1; + break; + } + h1 = next(h1, tm1); + } + if (h_to_split != BGT::null_halfedge()) + { + halfedge_descriptor hnew = split_edge_and_incident_faces(h_to_split, get(vpm2, vd), tm1, vpm1); + put(matching_vertices_1, target(hnew, tm1), true); + put(matching_vertices_2, vd, true); + break; + } + break; // TODO: decide what to do here. snap on the closest face? + } + default: + { + // look for a vertex common to all faces + std::set vset; + halfedge_descriptor h1 = halfedge(closed_faces[0], tm1); + for (int i=0;i<3; ++i) + { + vset.insert(target(h1, tm1)); + h1 = next(h1, tm1); + } + for (auto fit = std::next(closed_faces.begin()); fit!=closed_faces.end(); ++fit) + { + h1 = halfedge(*fit, tm1); + std::set common; + for (int i=0;i<3; ++i) + { + if (vset.count(target(h1, tm1))==1) + common.insert(target(h1,tm1)); + h1 = next(h1, tm1); + } + common.swap(vset); + if (vset.empty()) break; + } + if (vset.size()!=1) break; // TODO snap on the closest? + //snap on the vertex + put(vpm2, vd, get(vpm1, *vset.begin())); + put(matching_vertices_1, *vset.begin(), true); + put(matching_vertices_2, vd, true); + } + } + } +} + +// TODO test on Thingi10K/raw_meshes/162100.stl after splitting it +template +void surface_snapping(TriangleMesh& tm1, TriangleMesh& tm2, double epsilon) +{ + using VPM = typename property_map_selector::type; + using K = typename Kernel_traits::value_type>::type; + + using AABB_primitive_ = AABB_face_graph_triangle_primitive; + using AABB_traits_ = AABB_traits; + using AABB_tree_ = AABB_tree; + + using BGT = boost::graph_traits; + using vertex_descriptor = typename BGT::vertex_descriptor; + + VPM vpm1 = get(boost::vertex_point, tm1); + VPM vpm2 = get(boost::vertex_point, tm2); + +// check for identical points in tm1 and tm2 + // Note does not work well in parallel: bool --> int + auto matching_vertices_1 = get(CGAL::dynamic_vertex_property_t(), tm1); + auto matching_vertices_2 = get(CGAL::dynamic_vertex_property_t(), tm2); + std::map pts_map; + + for (vertex_descriptor vd : vertices(tm1)) + { + auto insert_res = pts_map.emplace(get(vpm1, vd), vd); + if (!insert_res.second) + std::cerr << "Warning duplicated vertex in tm1, better solve that first!\n"; + } + + for (vertex_descriptor vd : vertices(tm2)) + { + auto insert_res = pts_map.emplace(get(vpm2, vd), BGT::null_vertex()); + if (!insert_res.second) + { + if (insert_res.first->second == BGT::null_vertex()) + std::cerr << "Warning duplicated vertex in tm2, better solve that first!\n"; + else + { + if (get(matching_vertices_1, insert_res.first->second)) + { + std::cerr << "Two vertices of tm2 closed to the same vertex of tm1\n"; + } + put(matching_vertices_1, insert_res.first->second, true); + put(matching_vertices_2, vd, true); + } + } + } + + AABB_tree_ tm1_tree(faces(tm1).begin(), faces(tm1).end(), tm1); + AABB_tree_ tm2_tree(faces(tm2).begin(), faces(tm2).end(), tm2); + + double eps2 = epsilon * epsilon; + surface_snapping_impl_tm1_on_tm2(tm1_tree, tm1, tm2, vpm1, vpm2, + matching_vertices_1, matching_vertices_2, + eps2); + surface_snapping_impl_tm1_on_tm2(tm2_tree, tm2, tm1, vpm2, vpm1, + matching_vertices_2, matching_vertices_1, + eps2); + + + + // TODO: apply the same thing using tm2_tree ( + +} + +} } } /// end of CGAL::Polygon_mesh_processing::experimental + +#endif // CGAL_POLYGON_MESH_PROCESSING_INTERNAL_SURFACE_SNAPPING_H \ No newline at end of file diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt index 44e32386664c..b9e21bf8fab5 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt @@ -37,6 +37,7 @@ create_single_source_cgal_program("test_corefinement_and_constraints.cpp") create_single_source_cgal_program("test_corefinement_and_constraints_nearest.cpp") create_single_source_cgal_program("test_corefinement_bool_op.cpp") create_single_source_cgal_program("test_corefine.cpp") +create_single_source_cgal_program("snap_and_corefinement_union.cpp") create_single_source_cgal_program("test_coref_epic_points_identity.cpp") create_single_source_cgal_program("test_does_bound_a_volume.cpp") create_single_source_cgal_program("test_pmp_clip.cpp") diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/snap_and_corefinement_union.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/snap_and_corefinement_union.cpp new file mode 100644 index 000000000000..69192938224c --- /dev/null +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/snap_and_corefinement_union.cpp @@ -0,0 +1,44 @@ +#include +#include + +#include +#include +#include + +#include +#include + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Surface_mesh Mesh; + +namespace PMP = CGAL::Polygon_mesh_processing; + +int main(int argc, char* argv[]) +{ + const std::string filename1 = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/blobby.off"); + const std::string filename2 = (argc > 2) ? argv[2] : CGAL::data_file_path("meshes/eight.off"); + const double epsilon = (argc > 3) ? atof(argv[3]) : 0.1; + + Mesh mesh1, mesh2; + if(!PMP::IO::read_polygon_mesh(filename1, mesh1) || !PMP::IO::read_polygon_mesh(filename2, mesh2)) + { + std::cerr << "Invalid input." << std::endl; + return 1; + } + + PMP::experimental::surface_snapping(mesh1, mesh2, epsilon); + + Mesh out; + bool valid_union = PMP::corefine_and_compute_union(mesh1,mesh2, out); + + if(valid_union) + { + std::cout << "Union was successfully computed\n"; + CGAL::IO::write_polygon_mesh("snap_union.off", out, CGAL::parameters::stream_precision(17)); + return 0; + } + + std::cout << "Union could not be computed\n"; + + return 1; +}