From cdc3bd22cf6218ecfdee1ea8d092f7afa86bfbab Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Tue, 29 Nov 2022 12:31:04 +0100 Subject: [PATCH] Fix project_on_surface's projection direction See also https://github.com/CGAL/cgal/pull/5209 --- Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h | 24 ++++----- Mesh_3/include/CGAL/Mesh_triangulation_3.h | 6 ++- .../CGAL/Periodic_3_mesh_triangulation_3.h | 50 ++++++------------- 3 files changed, 28 insertions(+), 52 deletions(-) diff --git a/Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h b/Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h index 6a5e0923a01..2844dfd5a63 100644 --- a/Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h +++ b/Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h @@ -3452,17 +3452,10 @@ get_least_square_surface_plane(const Vertex_handle& v, (patch_index == Surface_patch_index() || c3t3_.surface_patch_index(f) == patch_index) ) { - ref_facet = f; + if(ref_facet.first == Cell_handle()) + ref_facet = f; - // In the case of a periodic triangulation, the incident facets of a point - // do not necessarily have the same offsets. Worse, the surface centers - // might not have the same offset as their facet. Thus, no solution except - // calling a function 'get_closest_triangle(p, t)' that simply returns t - // for a non-periodic triangulation, and checks all possible offsets for - // periodic triangulations - - Triangle t = c3t3_.triangulation().triangle(f); - Triangle ct = tr_.get_closest_triangle(cp(position), t); + const Triangle& ct = tr_.get_incident_triangle(f, v); triangles.push_back(ct); } } @@ -3474,7 +3467,6 @@ get_least_square_surface_plane(const Vertex_handle& v, // Compute least square fitting plane Plane_3 plane; Bare_point point; - CGAL::linear_least_squares_fitting_3(triangles.begin(), triangles.end(), plane, @@ -3483,12 +3475,13 @@ get_least_square_surface_plane(const Vertex_handle& v, tr_.geom_traits(), Default_diagonalize_traits()); - return std::make_pair(plane, - ref_facet.first->get_facet_surface_center(ref_facet.second)); + + // The surface center of a facet might have an offset in periodic triangulations + const Bare_point& ref_facet_scp = ref_facet.first->get_facet_surface_center(ref_facet.second); + const Bare_point& ref_point = tr_.get_closest_point(cp(position), ref_facet_scp); + return std::make_pair(plane, ref_point); } - - template typename C3T3_helpers::Bare_point C3T3_helpers:: @@ -3510,6 +3503,7 @@ project_on_surface_if_possible(const Vertex_handle& v, const Bare_point& p, Surface_patch_index index) const { + // @todo should call below if it's available... // return domain_.project_on_surface(p); typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object(); diff --git a/Mesh_3/include/CGAL/Mesh_triangulation_3.h b/Mesh_3/include/CGAL/Mesh_triangulation_3.h index 8d8012806c7..2970dc92b58 100644 --- a/Mesh_3/include/CGAL/Mesh_triangulation_3.h +++ b/Mesh_3/include/CGAL/Mesh_triangulation_3.h @@ -62,12 +62,14 @@ public: typedef typename Base::Triangle Triangle; typedef typename Base::Vertex_handle Vertex_handle; + typedef typename Base::Facet Facet; typedef typename Base::Cell_handle Cell_handle; typedef typename Geom_traits::Vector_3 Vector; using Base::geom_traits; using Base::point; + using Base::triangle; static std::string io_signature() { return Get_io_signature()(); } @@ -83,9 +85,9 @@ public: return q; } - const Triangle& get_closest_triangle(const Bare_point& /*p*/, const Triangle& t) const + Triangle get_incident_triangle(const Facet& f, const Vertex_handle) const { - return t; + return triangle(f); } void set_point(const Vertex_handle v, diff --git a/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h b/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h index 7b3eefe9105..a9a952215ec 100644 --- a/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h +++ b/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h @@ -468,45 +468,25 @@ public: return cwp(get_closest_point(cp(wp), cp(wq)), cw(wq)); } - Triangle get_closest_triangle(const Bare_point& p, const Triangle& t) const + // returns the triangle corresponding to f, with a geometric shift + // so that it is incident to ref_v's canonical position + Triangle get_incident_triangle(const Facet& f, const Vertex_handle ref_v) const { - typename Geom_traits::Construct_vector_3 cv = geom_traits().construct_vector_3_object(); + typename Geom_traits::Construct_point_3 cp = geom_traits().construct_point_3_object(); typename Geom_traits::Construct_translated_point_3 tr = geom_traits().construct_translated_point_3_object(); - typename Geom_traits::Compute_squared_distance_3 csd = geom_traits().compute_squared_distance_3_object(); + typename Geom_traits::Construct_vector_3 cv = geom_traits().construct_vector_3_object(); + typename Geom_traits::Construct_triangle_3 ct = geom_traits().construct_triangle_3_object(); - // It doesn't matter which point we use to canonicalize the triangle since we have to look - // at all the neighboring copies anyway since we do not have control of 'p'. - Bare_point canon_p0 = canonicalize_point(t[0]); - Vector_3 move_to_canonical = cv(t[0], canon_p0); - const std::array ct = { canon_p0, - tr(t[1], move_to_canonical), - tr(t[2], move_to_canonical) }; + CGAL_precondition(f.first->has_vertex(ref_v)); + const int ref_v_pos = f.first->index(ref_v); + const Bare_point& ref_p = cp(point(ref_v)); + const Bare_point& ref_p_in_f = cp(point(f.first, ref_v_pos)); + Vector_3 move_to_canonical = cv(ref_p_in_f, ref_p); - FT min_sq_dist = std::numeric_limits::infinity(); - Triangle rt; - for(int i = -1; i < 2; ++i) { - for(int j = -1; j < 2; ++j) { - for(int k = -1; k < 2; ++k) { - - const Triangle tt( - construct_point(std::make_pair(ct[0], Offset(i, j, k))), - construct_point(std::make_pair(ct[1], Offset(i, j, k))), - construct_point(std::make_pair(ct[2], Offset(i, j, k)))); - - const FT sq_dist = csd(p, tt); - - if(sq_dist == FT(0)) - return rt; - - if(sq_dist < min_sq_dist) { - rt = tt; - min_sq_dist = sq_dist; - } - } - } - } - - return rt; + const int s = f.second; + return ct(tr(cp(point(f.first, (s+1)%4)), move_to_canonical), + tr(cp(point(f.first, (s+2)%4)), move_to_canonical), + tr(cp(point(f.first, (s+3)%4)), move_to_canonical)); } // Warning: This is a periodic version that computes the smallest possible