From 0c8739a16ed1a82b1b3239eb3ce9ec91b41b4106 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Mon, 8 Jul 2019 15:39:33 +0200 Subject: [PATCH 01/39] Re-order convenience overloads to avoid forward declarations in upcoming code --- .../Polygon_mesh_processing/compute_normal.h | 108 +++++++++--------- 1 file changed, 57 insertions(+), 51 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h index a0039b0273b..87e4c73fc27 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h @@ -49,7 +49,7 @@ namespace internal { { typename GT::FT norm = CGAL::approximate_sqrt( traits.compute_squared_length_3_object()(v)); - //If the vector is small enough, approx_sqrt might return 0, and then we get nan values. + //If the vector is small enough, approx_sqrt might return 0, and then we get nan values. //To avoid that, we check the resulted norm. If it is 0, we don't normalize. if(norm != 0) { @@ -67,7 +67,7 @@ namespace internal { traits.construct_vector_3_object()(p1, p2), traits.construct_vector_3_object()(p1, p0)); - //cross-product(AB, AC)'s norm is the area of the parallelogram + //cross-product(AB, AC)'s norm is the area of the parallelogram //formed by these 2 vectors. //the triangle's area is half of it return traits.construct_scaled_vector_3_object()(n, 0.5); @@ -159,6 +159,20 @@ compute_face_normal(typename boost::graph_traits::face_descriptor f return normal; } +///\cond SKIP_IN_MANUAL + +// compute_face_normal overload +template +typename GetGeomTraits::type::Vector_3 +compute_face_normal(typename boost::graph_traits::face_descriptor f, + const PolygonMesh& pmesh) +{ + return compute_face_normal(f, pmesh, + CGAL::Polygon_mesh_processing::parameters::all_default()); +} + +/// \endcond + /** * \ingroup PMP_normal_grp * computes the outward unit vector normal for all faces of the polygon mesh. @@ -198,6 +212,17 @@ compute_face_normals(const PolygonMesh& pmesh } } +///\cond SKIP_IN_MANUAL + +// compute_face_normals overload +template +void compute_face_normals(const PolygonMesh& pmesh, FaceNormalMap fnm) +{ + compute_face_normals(pmesh, fnm, CGAL::Polygon_mesh_processing::parameters::all_default()); +} + +/// \endcond + /** * \ingroup PMP_normal_grp * computes the unit normal at vertex `v` as the average of the normals of incident faces. @@ -275,6 +300,19 @@ compute_vertex_normal(typename boost::graph_traits::vertex_descript return normal; } +///\cond SKIP_IN_MANUAL + +// compute_vertex_normal overloads +template +typename GetGeomTraits::type::Vector_3 +compute_vertex_normal(typename boost::graph_traits::vertex_descriptor v, + const PolygonMesh& pmesh) +{ + return compute_vertex_normal(v, pmesh, CGAL::Polygon_mesh_processing::parameters::all_default()); +} + +/// \endcond + /** * \ingroup PMP_normal_grp * computes the outward unit vector normal for all vertices of the polygon mesh. @@ -316,6 +354,17 @@ compute_vertex_normals(const PolygonMesh& pmesh } } +///\cond SKIP_IN_MANUAL + +// compute_vertex_normals overloads +template +void compute_vertex_normals(const PolygonMesh& pmesh, VertexNormalMap vnm) +{ + compute_vertex_normals(pmesh, vnm, CGAL::Polygon_mesh_processing::parameters::all_default()); +} + +/// \endcond + /** * \ingroup PMP_normal_grp * computes the outward unit vector normal for all vertices and faces of the polygon mesh. @@ -360,62 +409,19 @@ compute_normals(const PolygonMesh& pmesh } ///\cond SKIP_IN_MANUAL -// compute_vertex_normal overloads -template -typename CGAL::Kernel_traits< typename property_map_value::type>::Kernel::Vector_3 -compute_vertex_normal( - typename boost::graph_traits::vertex_descriptor v, - const PolygonMesh& pmesh) -{ - return compute_vertex_normal(v, pmesh, - CGAL::Polygon_mesh_processing::parameters::all_default()); -} - -// compute_vertex_normals overloads -template -void -compute_vertex_normals(const PolygonMesh& pmesh, - VertexNormalMap vnm) -{ - compute_vertex_normals(pmesh, vnm, - CGAL::Polygon_mesh_processing::parameters::all_default()); -} - -// compute_face_normal overload -template -typename CGAL::Kernel_traits < typename property_map_value::type>::Kernel::Vector_3 -compute_face_normal( - typename boost::graph_traits::face_descriptor f, - const PolygonMesh& pmesh) -{ - return compute_face_normal(f, pmesh, - CGAL::Polygon_mesh_processing::parameters::all_default()); -} - -// compute_face_normals overload -template -void -compute_face_normals(const PolygonMesh& pmesh, FaceNormalMap fnm) -{ - compute_face_normals(pmesh, fnm, - CGAL::Polygon_mesh_processing::parameters::all_default()); -} // compute_normals overload template -void -compute_normals(const PolygonMesh& pmesh, - VertexNormalMap vnm, - FaceNormalMap fnm) +void compute_normals(const PolygonMesh& pmesh, + VertexNormalMap vnm, + FaceNormalMap fnm) { - compute_normals(pmesh, vnm, fnm, - CGAL::Polygon_mesh_processing::parameters::all_default()); + compute_normals(pmesh, vnm, fnm, CGAL::Polygon_mesh_processing::parameters::all_default()); } /// \endcond -} - -} // end of namespace CGAL::Polygon_mesh_processing +} // namespace Polygon_mesh_processing +} // namespace CGAL #endif // CGAL_POLYGON_MESH_PROCESSING_COMPUTE_NORMAL_H From 401aea2500da826facdea61ac0c2bc0890dbd0ac Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Mon, 8 Jul 2019 15:40:08 +0200 Subject: [PATCH 02/39] Add new method to compute vertex normals --- .../Polygon_mesh_processing/compute_normal.h | 516 ++++++++++++++++++ 1 file changed, 516 insertions(+) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h index 87e4c73fc27..5bb1733778a 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h @@ -223,6 +223,522 @@ void compute_face_normals(const PolygonMesh& pmesh, FaceNormalMap fnm) /// \endcond + +// @TMP +// ----------------------------------------------------------------- +template +typename K::FT fix_sine(typename K::FT sine) +{ + if(sine >= 1) + return 1; + else if(sine <= -1) + return -1; + else + return sine; +} + +template +typename K::FT compute_angle_rad(const typename K::Vector_3& u, + const typename K::Vector_3& v) +{ + typedef typename K::FT NT; + + // check + NT product = CGAL::approximate_sqrt(u * u) * CGAL::approximate_sqrt(v * v); + if(product == 0) + return 0; + + // cosine + NT dot = (u * v); + NT cosine = dot / product; + + return std::acos(fix_sine(cosine)); +} + +// -> -> +// Returns the angle (in radians) of (P,Q,R) corner (i.e. QP, QR angle). +template +typename K::FT compute_angle_rad(const typename K::Point_3& P, + const typename K::Point_3& Q, + const typename K::Point_3& R) +{ + typedef typename K::Vector_3 Vector_3; + + Vector_3 u = P - Q; + Vector_3 v = R - Q; + + return compute_angle_rad(u, v); +} +// ----------------------------------------------------------------- + +template +bool almost_equal(const typename GT::Vector_3 v1, const typename GT::Vector_3 v2) +{ + return (CGAL::abs(1 - GT().compute_scalar_product_3_object()(v1, v2)) < 1e-10); +} + +template +bool +does_enclose_other_normals(const int i, const int j, const int k, + const typename K::Vector_3& nb, + const typename K::FT sp_bi, + std::vector::face_descriptor> incident_faces, + const FaceNormalVector& normalized_face_normals, + const K& traits) +{ + typedef typename K::FT FT; + typedef typename K::Vector_3 Vector_3; + + typename K::Compute_scalar_product_3 sp = traits.compute_scalar_product_3_object(); + + // check that this min circle defined by the diameter contains the other points + const std::size_t nif = incident_faces.size(); + for(int l=0; l +typename GT::Vector_3 compute_normals_bisector(const typename GT::Vector_3 ni, + const typename GT::Vector_3 nj, + const GT& traits) +{ + if(ni == nj) + return ni; + + typename GT::Vector_3 nb = traits.construct_sum_of_vectors_3_object()(ni, nj); + CGAL::Polygon_mesh_processing::internal::normalize(nb, traits); // prob not necessary @fixme + + return nb; +} + +template +std::pair +compute_most_visible_normal_2_points(std::vector::face_descriptor> incident_faces, + const FaceNormalVector& normalized_face_normals, + const GT& traits) +{ + typedef typename GT::FT FT; + typedef typename GT::Vector_3 Vector_3; + + typename GT::Compute_scalar_product_3 sp = traits.compute_scalar_product_3_object(); + + FT min_sp = -1; + + Vector_3 n = CGAL::NULL_VECTOR; + + const std::size_t nif = incident_faces.size(); + for(int i=0; i= 0); + if(sp_bi <= min_sp) + continue; + + if(!does_enclose_other_normals(i, j, -1 /*NA*/, nb, sp_bi, incident_faces, normalized_face_normals, traits)) + continue; + + min_sp = sp_bi; + n = nb; + } + } + + return std::make_pair(n, (n != CGAL::NULL_VECTOR)); +} + +template +typename GT::Vector_3 +compute_most_visible_normal_3_points(std::vector::face_descriptor> incident_faces, + const FaceNormalVector& normalized_face_normals, + const GT& traits) +{ + typedef typename GT::FT FT; + typedef typename GT::Point_3 Point_3; + typedef typename GT::Vector_3 Vector_3; + + typename GT::Compute_scalar_product_3 sp = traits.compute_scalar_product_3_object(); + typename GT::Construct_sum_of_vectors_3 csv = traits.construct_sum_of_vectors_3_object(); + + FT min_sp = -1; + + Vector_3 n = CGAL::NULL_VECTOR; + + const std::size_t nif = incident_faces.size(); + + for(int i=0; i(ni, nj)) + { + if(almost_equal(nj, nk)) + nb = ni; + else // ni == nj, but nj != nk + nb = compute_normals_bisector(nj, nk, traits); + } + else if(almost_equal(ni, nk)) // ni != nj + { + nb = compute_normals_bisector(nj, nk, traits); + } + else if(almost_equal(nj, nk)) // ni != nj + { + nb = compute_normals_bisector(ni, nk, traits); + } + else + { + CGAL_assertion(ni != nj); + CGAL_assertion(ni != nk); + CGAL_assertion(nj != nk); + CGAL_assertion(!CGAL::collinear(CGAL::ORIGIN + ni, CGAL::ORIGIN + nj, CGAL::ORIGIN + nk)); + + const Point_3 c = CGAL::circumcenter(CGAL::ORIGIN + ni, CGAL::ORIGIN + nj, CGAL::ORIGIN + nk); + CGAL_assertion(c != CGAL::ORIGIN); + + nb = traits.construct_vector_3_object()(CGAL::ORIGIN, c); + CGAL::Polygon_mesh_processing::internal::normalize(nb, traits); // prob not necessary @fixme + } + + CGAL_assertion(nb != CGAL::NULL_VECTOR); + FT sp_bi = sp(nb, ni); + + if(sp_bi < 0) + { + nb = traits.construct_opposite_vector_3_object()(nb); + sp_bi = - sp_bi; + } + + if(sp_bi <= min_sp) + continue; + + if(!does_enclose_other_normals(i, j, k, nb, sp_bi, incident_faces, normalized_face_normals, traits)) + continue; + + min_sp = sp_bi; + n = nb; + } + } + } + + CGAL_assertion(n != CGAL::NULL_VECTOR); + return n; +} + +// Complexity is high, but valence is usually low, so that's ok +template +typename GT::Vector_3 +compute_vertex_normal_most_visible_min_circle(typename boost::graph_traits::vertex_descriptor vd, + const FaceNormalVector& normalized_face_normals, + const PolygonMesh& pmesh, + const GT& traits) +{ + typedef typename GT::FT FT; + typedef typename GT::Vector_3 Vector_3; + + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + typedef typename boost::graph_traits::face_descriptor face_descriptor; + +#ifdef CGAL_MOST_VISIBLE_NORMAL_VERBOSE + std::cout << std::endl << std::endl; + std::cout << "----------------------------------------------------------------------" << std::endl; + std::cout << "compute_vertex_normal_most_visible_min_circle at " << pmesh.point(vd) << std::endl; +#endif + + halfedge_descriptor hd = halfedge(vd, pmesh); + + std::vector incident_faces; + for(face_descriptor fd : CGAL::faces_around_target(hd, pmesh)) + { + if(fd == boost::graph_traits::null_face()) + continue; + + incident_faces.push_back(fd); + } + + std::pair res = compute_most_visible_normal_2_points(incident_faces, normalized_face_normals, traits); + if(res.second) + return res.first; + + return compute_most_visible_normal_3_points(incident_faces, normalized_face_normals, traits); +} + +template +typename GT::Vector_3 +compute_vertex_normal_most_visible_with_optimization(typename boost::graph_traits::vertex_descriptor vd, + const FaceNormalVector& normalized_face_normals, + const PolygonMesh& pmesh, + const GT& traits) +{ + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + typedef typename boost::graph_traits::face_descriptor face_descriptor; + + typedef typename GT::FT FT; + typedef typename GT::Vector_3 Vector_3; + + typename GT::Construct_scaled_vector_3 csclv = traits.construct_scaled_vector_3_object(); + typename GT::Construct_sum_of_vectors_3 csv = traits.construct_sum_of_vectors_3_object(); + typename GT::Compute_scalar_product_3 sp = traits.compute_scalar_product_3_object(); + +#ifdef CGAL_MOST_VISIBLE_NORMAL_VERBOSE + std::cout << "----------------------------------------------------------------------" << std::endl; + std::cout << "compute_vertex_normal_most_visible_with_optimization at " << pmesh.point(vd) << std::endl; +#endif + + halfedge_descriptor hd = halfedge(vd, pmesh); + + std::map weights; + Vector_3 normal = CGAL::NULL_VECTOR; + +// initial weight based on the number of incident faces +#ifdef CGAL_PMP_MOST_VISIBLE_NORMAL_INITIALIZE_WITH_EQUAL_WEIGHTS + typedef typename CGAL::Halfedge_around_target_iterator::difference_type difference_type; + difference_type incident_faces_n = 0; + + for(face_descriptor fd : CGAL::faces_around_target(hd, pmesh)) + { + if(fd == boost::graph_traits::null_face()) + continue; + + ++incident_faces_n; + } + + CGAL_assertion(incident_faces_n > 0); + + // Initial normal is just weighted by 1/n + const FT initial_weight = FT(1)/incident_faces_n; + + for(face_descriptor fd : CGAL::faces_around_target(hd, pmesh)) + { + if(fd == boost::graph_traits::null_face()) + continue; + + normal = csv(normal, csclv(normalized_face_normals.at(fd), initial_weight)); + weights[face] = initial_weight; + } +#else // angle-based initial weight + for(face_descriptor fd : CGAL::faces_around_target(hd, pmesh)) + { + if(fd == boost::graph_traits::null_face()) + continue; + + halfedge_descriptor can_hd = halfedge(fd, pmesh); + while(target(can_hd, pmesh) != vd) + can_hd = next(can_hd, pmesh); + + FT weight = compute_angle_rad(pmesh.point(target(next(can_hd, pmesh), pmesh)), + pmesh.point(target(can_hd, pmesh)), + pmesh.point(source(can_hd, pmesh))); + +// #define CGAL_MOST_VISIBLE_NORMAL_VERBOSE +#ifdef CGAL_MOST_VISIBLE_NORMAL_VERBOSE + std::cout << " weight: " << weight << std::endl; +#endif + + normal = csv(normal, csclv(normalized_face_normals.at(fd), weight)); + weights[fd] = weight; + } +#endif + + CGAL::Polygon_mesh_processing::internal::normalize(normal, traits); + + // Iterative process to find the most visible + bool converged = false; + int iter = 0; + FT old_sp = std::numeric_limits::infinity(); + for(;;) + { +#ifdef CGAL_MOST_VISIBLE_NORMAL_VERBOSE + std::cout << "Current best normal: " << normal << " for vertex: " << pmesh.point(vd) << std::endl; +#endif + + ++iter; + + // compute scaling coefficients + boost::container::flat_map alphas; + alphas.reserve(16); + FT alpha_sum = 0; + + for(face_descriptor fd : CGAL::faces_around_target(hd, pmesh)) + { + if(fd == boost::graph_traits::null_face()) + continue; + + FT normals_sp = sp(normal, normalized_face_normals.at(fd)); + if(normals_sp > 1) // everything is normalized so that shouldn't happen in theory, but it does in practice + normals_sp = 1; + if(normals_sp < -1) + normals_sp = -1; + + const FT alpha = std::acos(normals_sp); + +#ifdef CGAL_MOST_VISIBLE_NORMAL_VERBOSE + std::cout << " alpha: " << alpha << " (" << normalized_face_normals.at(fd) << ")" << std::endl; +#endif + + CGAL_assertion(alpha >= 0); + alphas[fd] = alpha; + alpha_sum += alpha; + } + + if(alpha_sum == 0) // all the sp are 1 and thus all the vectors are identical + return normal; + + CGAL_assertion(weights.size() == alphas.size()); + + // compute the new weights + FT weighted_sum = 0; + + for(face_descriptor fd : CGAL::faces_around_target(hd, pmesh)) + { + if(fd == boost::graph_traits::null_face()) + continue; + + FT& weight = weights[fd]; + weight *= alphas.at(fd) / alpha_sum; + +#ifdef CGAL_MOST_VISIBLE_NORMAL_VERBOSE + std::cout << " weights[" << fd << "] = " << weight << std::endl; +#endif + + weighted_sum += weight; + } + + if(weighted_sum == 0) // not too sure about that one... @fixme (correct if all weights are 0, but what otherwise?) + return normal; + + for(face_descriptor fd : CGAL::faces_around_target(hd, pmesh)) + { + if(fd == boost::graph_traits::null_face()) + continue; + + weights[fd] /= weighted_sum; + } + + // compute the new normal + Vector_3 new_normal_base = CGAL::NULL_VECTOR; + for(face_descriptor fd : CGAL::faces_around_target(hd, pmesh)) + { + if(fd == boost::graph_traits::null_face()) + continue; + +#ifdef CGAL_MOST_VISIBLE_NORMAL_VERBOSE + std::cout << "add: " << normalized_face_normals.at(fd) << " with weight: " << weights[fd] << std::endl; +#endif + new_normal_base = csv(new_normal_base, csclv(normalized_face_normals.at(fd), weights[fd])); + } + + CGAL::Polygon_mesh_processing::internal::normalize(new_normal_base, traits); + + // check convergence + const FT bound = 1e-10; + const FT new_sp = sp(normal, new_normal_base); + converged = (CGAL::abs(new_sp - old_sp) < bound); + +#ifdef CGAL_MOST_VISIBLE_NORMAL_VERBOSE + std::cout << "new normal (base): " << new_normal_base << std::endl; + std::cout << " convergence ? " << old_sp << " " << new_sp << " diff: " << CGAL::abs(new_sp - old_sp) << std::endl; +#endif + + if(converged || iter > 1000000) // @fixme + { + normal = new_normal_base; + break; + } + + // a bit of relaxation + const FT beta = 0.01; + const Vector_3 beta_scaled_normal = csclv(new_normal_base, beta); + normal = csv(beta_scaled_normal, csclv(normal, 1 - beta)); + old_sp = new_sp; + } + +#ifdef CGAL_MOST_VISIBLE_NORMAL_VERBOSE + std::cout << "Converged in: " << iter << " iterations" << std::endl; +#endif + + return normal; +} + +template +void compute_most_visible_vertex_normals(VertexNormalMap& vertex_normal_map, + const PolygonMesh& pmesh, + const NamedParameters& np) +{ + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::face_descriptor face_descriptor; + + CGAL_precondition(CGAL::is_triangle_mesh(pmesh)); + + typedef typename GetVertexPointMap::const_type VPMap; +// VPMap vpmap = choose_param(get_param(np, internal_np::vertex_point), +// get_const_property_map(vertex_point, pmesh)); + + typedef typename GetGeomTraits::type GT; + const GT traits = choose_param(get_param(np, internal_np::geom_traits), GT()); + + typedef typename GT::Vector_3 Vector; + + std::cout << "compute face normals" << std::endl; + + std::vector normalized_face_normals(num_faces(pmesh)); + for(face_descriptor fd : faces(pmesh)) + { + Vector en = compute_face_normal(fd, pmesh); + CGAL::Polygon_mesh_processing::internal::normalize(en, traits); + normalized_face_normals[fd] = en; + } + + std::ofstream out("computed_normals.cgal.polylines.txt"); + + // we start by evaluating the translation normal for each vertex + for(vertex_descriptor vd : vertices(pmesh)) + { +#if 0 + Vector n = compute_vertex_normal_most_visible_with_optimization(vd, normalized_face_normals, pmesh, traits); +#else + Vector n = compute_vertex_normal_most_visible_min_circle(vd, normalized_face_normals, pmesh, traits); +#endif + + CGAL::Polygon_mesh_processing::internal::normalize(n, traits); + vertex_normal_map[vd] = n; + + out << "2 " << pmesh.point(vd) << " " << pmesh.point(vd) + 0.01 * n << std::endl; + } +} + +template +void compute_most_visible_vertex_normals(VertexNormalMap& vertex_normal_map, const PolygonMesh& pmesh) +{ + return compute_most_visible_vertex_normals(vertex_normal_map, pmesh, CGAL::parameters::all_default()); +} + /** * \ingroup PMP_normal_grp * computes the unit normal at vertex `v` as the average of the normals of incident faces. From eb086c7847c9ff568fe61dd03748b9d7899b5d38 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Tue, 9 Jul 2019 18:07:11 +0200 Subject: [PATCH 03/39] Fix border case of normal computation at a vertex with valence 1 --- .../include/CGAL/Polygon_mesh_processing/compute_normal.h | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h index 5bb1733778a..fdd87b54ee7 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h @@ -482,10 +482,15 @@ compute_vertex_normal_most_visible_min_circle(typename boost::graph_traits res = compute_most_visible_normal_2_points(incident_faces, normalized_face_normals, traits); if(res.second) return res.first; + CGAL_assertion(incident_faces.size() > 2); + return compute_most_visible_normal_3_points(incident_faces, normalized_face_normals, traits); } From f019ace14c47ac135cd411494bcf694651b965e0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Wed, 10 Jul 2019 08:57:21 +0200 Subject: [PATCH 04/39] Remove obsolete code --- .../Polygon_mesh_processing/compute_normal.h | 251 +----------------- 1 file changed, 1 insertion(+), 250 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h index fdd87b54ee7..13053a6d6c7 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h @@ -223,54 +223,6 @@ void compute_face_normals(const PolygonMesh& pmesh, FaceNormalMap fnm) /// \endcond - -// @TMP -// ----------------------------------------------------------------- -template -typename K::FT fix_sine(typename K::FT sine) -{ - if(sine >= 1) - return 1; - else if(sine <= -1) - return -1; - else - return sine; -} - -template -typename K::FT compute_angle_rad(const typename K::Vector_3& u, - const typename K::Vector_3& v) -{ - typedef typename K::FT NT; - - // check - NT product = CGAL::approximate_sqrt(u * u) * CGAL::approximate_sqrt(v * v); - if(product == 0) - return 0; - - // cosine - NT dot = (u * v); - NT cosine = dot / product; - - return std::acos(fix_sine(cosine)); -} - -// -> -> -// Returns the angle (in radians) of (P,Q,R) corner (i.e. QP, QR angle). -template -typename K::FT compute_angle_rad(const typename K::Point_3& P, - const typename K::Point_3& Q, - const typename K::Point_3& R) -{ - typedef typename K::Vector_3 Vector_3; - - Vector_3 u = P - Q; - Vector_3 v = R - Q; - - return compute_angle_rad(u, v); -} -// ----------------------------------------------------------------- - template bool almost_equal(const typename GT::Vector_3 v1, const typename GT::Vector_3 v2) { @@ -494,203 +446,6 @@ compute_vertex_normal_most_visible_min_circle(typename boost::graph_traits(incident_faces, normalized_face_normals, traits); } -template -typename GT::Vector_3 -compute_vertex_normal_most_visible_with_optimization(typename boost::graph_traits::vertex_descriptor vd, - const FaceNormalVector& normalized_face_normals, - const PolygonMesh& pmesh, - const GT& traits) -{ - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - typedef typename boost::graph_traits::face_descriptor face_descriptor; - - typedef typename GT::FT FT; - typedef typename GT::Vector_3 Vector_3; - - typename GT::Construct_scaled_vector_3 csclv = traits.construct_scaled_vector_3_object(); - typename GT::Construct_sum_of_vectors_3 csv = traits.construct_sum_of_vectors_3_object(); - typename GT::Compute_scalar_product_3 sp = traits.compute_scalar_product_3_object(); - -#ifdef CGAL_MOST_VISIBLE_NORMAL_VERBOSE - std::cout << "----------------------------------------------------------------------" << std::endl; - std::cout << "compute_vertex_normal_most_visible_with_optimization at " << pmesh.point(vd) << std::endl; -#endif - - halfedge_descriptor hd = halfedge(vd, pmesh); - - std::map weights; - Vector_3 normal = CGAL::NULL_VECTOR; - -// initial weight based on the number of incident faces -#ifdef CGAL_PMP_MOST_VISIBLE_NORMAL_INITIALIZE_WITH_EQUAL_WEIGHTS - typedef typename CGAL::Halfedge_around_target_iterator::difference_type difference_type; - difference_type incident_faces_n = 0; - - for(face_descriptor fd : CGAL::faces_around_target(hd, pmesh)) - { - if(fd == boost::graph_traits::null_face()) - continue; - - ++incident_faces_n; - } - - CGAL_assertion(incident_faces_n > 0); - - // Initial normal is just weighted by 1/n - const FT initial_weight = FT(1)/incident_faces_n; - - for(face_descriptor fd : CGAL::faces_around_target(hd, pmesh)) - { - if(fd == boost::graph_traits::null_face()) - continue; - - normal = csv(normal, csclv(normalized_face_normals.at(fd), initial_weight)); - weights[face] = initial_weight; - } -#else // angle-based initial weight - for(face_descriptor fd : CGAL::faces_around_target(hd, pmesh)) - { - if(fd == boost::graph_traits::null_face()) - continue; - - halfedge_descriptor can_hd = halfedge(fd, pmesh); - while(target(can_hd, pmesh) != vd) - can_hd = next(can_hd, pmesh); - - FT weight = compute_angle_rad(pmesh.point(target(next(can_hd, pmesh), pmesh)), - pmesh.point(target(can_hd, pmesh)), - pmesh.point(source(can_hd, pmesh))); - -// #define CGAL_MOST_VISIBLE_NORMAL_VERBOSE -#ifdef CGAL_MOST_VISIBLE_NORMAL_VERBOSE - std::cout << " weight: " << weight << std::endl; -#endif - - normal = csv(normal, csclv(normalized_face_normals.at(fd), weight)); - weights[fd] = weight; - } -#endif - - CGAL::Polygon_mesh_processing::internal::normalize(normal, traits); - - // Iterative process to find the most visible - bool converged = false; - int iter = 0; - FT old_sp = std::numeric_limits::infinity(); - for(;;) - { -#ifdef CGAL_MOST_VISIBLE_NORMAL_VERBOSE - std::cout << "Current best normal: " << normal << " for vertex: " << pmesh.point(vd) << std::endl; -#endif - - ++iter; - - // compute scaling coefficients - boost::container::flat_map alphas; - alphas.reserve(16); - FT alpha_sum = 0; - - for(face_descriptor fd : CGAL::faces_around_target(hd, pmesh)) - { - if(fd == boost::graph_traits::null_face()) - continue; - - FT normals_sp = sp(normal, normalized_face_normals.at(fd)); - if(normals_sp > 1) // everything is normalized so that shouldn't happen in theory, but it does in practice - normals_sp = 1; - if(normals_sp < -1) - normals_sp = -1; - - const FT alpha = std::acos(normals_sp); - -#ifdef CGAL_MOST_VISIBLE_NORMAL_VERBOSE - std::cout << " alpha: " << alpha << " (" << normalized_face_normals.at(fd) << ")" << std::endl; -#endif - - CGAL_assertion(alpha >= 0); - alphas[fd] = alpha; - alpha_sum += alpha; - } - - if(alpha_sum == 0) // all the sp are 1 and thus all the vectors are identical - return normal; - - CGAL_assertion(weights.size() == alphas.size()); - - // compute the new weights - FT weighted_sum = 0; - - for(face_descriptor fd : CGAL::faces_around_target(hd, pmesh)) - { - if(fd == boost::graph_traits::null_face()) - continue; - - FT& weight = weights[fd]; - weight *= alphas.at(fd) / alpha_sum; - -#ifdef CGAL_MOST_VISIBLE_NORMAL_VERBOSE - std::cout << " weights[" << fd << "] = " << weight << std::endl; -#endif - - weighted_sum += weight; - } - - if(weighted_sum == 0) // not too sure about that one... @fixme (correct if all weights are 0, but what otherwise?) - return normal; - - for(face_descriptor fd : CGAL::faces_around_target(hd, pmesh)) - { - if(fd == boost::graph_traits::null_face()) - continue; - - weights[fd] /= weighted_sum; - } - - // compute the new normal - Vector_3 new_normal_base = CGAL::NULL_VECTOR; - for(face_descriptor fd : CGAL::faces_around_target(hd, pmesh)) - { - if(fd == boost::graph_traits::null_face()) - continue; - -#ifdef CGAL_MOST_VISIBLE_NORMAL_VERBOSE - std::cout << "add: " << normalized_face_normals.at(fd) << " with weight: " << weights[fd] << std::endl; -#endif - new_normal_base = csv(new_normal_base, csclv(normalized_face_normals.at(fd), weights[fd])); - } - - CGAL::Polygon_mesh_processing::internal::normalize(new_normal_base, traits); - - // check convergence - const FT bound = 1e-10; - const FT new_sp = sp(normal, new_normal_base); - converged = (CGAL::abs(new_sp - old_sp) < bound); - -#ifdef CGAL_MOST_VISIBLE_NORMAL_VERBOSE - std::cout << "new normal (base): " << new_normal_base << std::endl; - std::cout << " convergence ? " << old_sp << " " << new_sp << " diff: " << CGAL::abs(new_sp - old_sp) << std::endl; -#endif - - if(converged || iter > 1000000) // @fixme - { - normal = new_normal_base; - break; - } - - // a bit of relaxation - const FT beta = 0.01; - const Vector_3 beta_scaled_normal = csclv(new_normal_base, beta); - normal = csv(beta_scaled_normal, csclv(normal, 1 - beta)); - old_sp = new_sp; - } - -#ifdef CGAL_MOST_VISIBLE_NORMAL_VERBOSE - std::cout << "Converged in: " << iter << " iterations" << std::endl; -#endif - - return normal; -} - template void compute_most_visible_vertex_normals(VertexNormalMap& vertex_normal_map, const PolygonMesh& pmesh, @@ -725,16 +480,12 @@ void compute_most_visible_vertex_normals(VertexNormalMap& vertex_normal_map, // we start by evaluating the translation normal for each vertex for(vertex_descriptor vd : vertices(pmesh)) { -#if 0 - Vector n = compute_vertex_normal_most_visible_with_optimization(vd, normalized_face_normals, pmesh, traits); -#else Vector n = compute_vertex_normal_most_visible_min_circle(vd, normalized_face_normals, pmesh, traits); -#endif CGAL::Polygon_mesh_processing::internal::normalize(n, traits); vertex_normal_map[vd] = n; - out << "2 " << pmesh.point(vd) << " " << pmesh.point(vd) + 0.01 * n << std::endl; + out << "2 " << pmesh.point(vd) << " " << pmesh.point(vd) + n << std::endl; } } From 9d1cec6be6736911c9df42dc35129fe23c4863f1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Thu, 18 Jul 2019 11:18:52 +0200 Subject: [PATCH 05/39] Fix indentation (cosmetic) --- .../Polygon_mesh_processing/compute_normal.h | 155 ++++++++---------- 1 file changed, 66 insertions(+), 89 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h index 13053a6d6c7..59e02ff731f 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h @@ -19,13 +19,11 @@ // // Author(s) : Jane Tournois - #ifndef CGAL_POLYGON_MESH_PROCESSING_COMPUTE_NORMAL_H #define CGAL_POLYGON_MESH_PROCESSING_COMPUTE_NORMAL_H #include - #include #include #include @@ -38,44 +36,39 @@ #include -namespace CGAL{ - -namespace Polygon_mesh_processing{ - +namespace CGAL { +namespace Polygon_mesh_processing { namespace internal { - template - void normalize(typename GT::Vector_3& v, const GT& traits) +template +void normalize(typename GT::Vector_3& v, const GT& traits) +{ + typename GT::FT norm = CGAL::approximate_sqrt(traits.compute_squared_length_3_object()(v)); + //If the vector is small enough, approx_sqrt might return 0, and then we get nan values. + //To avoid that, we check the resulted norm. If it is 0, we don't normalize. + if(norm != 0) { - typename GT::FT norm = CGAL::approximate_sqrt( - traits.compute_squared_length_3_object()(v)); - //If the vector is small enough, approx_sqrt might return 0, and then we get nan values. - //To avoid that, we check the resulted norm. If it is 0, we don't normalize. - if(norm != 0) - { - v = traits.construct_divided_vector_3_object()(v, norm ); - } - } - - template - typename GT::Vector_3 - triangle_normal(const Point& p0, const Point& p1, const Point& p2 - , const GT& traits) - { - typename GT::Vector_3 n = traits.construct_cross_product_vector_3_object()( - traits.construct_vector_3_object()(p1, p2), - traits.construct_vector_3_object()(p1, p0)); - - //cross-product(AB, AC)'s norm is the area of the parallelogram - //formed by these 2 vectors. - //the triangle's area is half of it - return traits.construct_scaled_vector_3_object()(n, 0.5); + v = traits.construct_divided_vector_3_object()(v, norm); } } -template +template +typename GT::Vector_3 +triangle_normal(const Point& p0, const Point& p1, const Point& p2, const GT& traits) +{ + typename GT::Vector_3 n = traits.construct_cross_product_vector_3_object()( + traits.construct_vector_3_object()(p1, p2), + traits.construct_vector_3_object()(p1, p0)); + + //cross-product(AB, AC)'s norm is the area of the parallelogram + //formed by these 2 vectors. + //the triangle's area is half of it + return traits.construct_scaled_vector_3_object()(n, 0.5); +} + +} // namespace internal + +template void sum_normals(const PM& pmesh, typename boost::graph_traits::face_descriptor f, VertexPointMap vpmap, @@ -88,9 +81,9 @@ void sum_normals(const PM& pmesh, halfedge_descriptor he = halfedge(f, pmesh); vertex_descriptor v = source(he, pmesh); const Point& pv = get(vpmap, v); - while (v != target(next(he, pmesh), pmesh)) + while(v != target(next(he, pmesh), pmesh)) { - const Point& pvn = get(vpmap, target(he, pmesh)); + const Point& pvn = get(vpmap, target(he, pmesh)); const Point& pvnn = get(vpmap, target(next(he, pmesh), pmesh)); Vector n = internal::triangle_normal(pv, pvn, pvnn, traits); @@ -133,9 +126,9 @@ Vector_3 #else typename GetGeomTraits::type::Vector_3 #endif -compute_face_normal(typename boost::graph_traits::face_descriptor f - , const PolygonMesh& pmesh - , const NamedParameters& np) +compute_face_normal(typename boost::graph_traits::face_descriptor f, + const PolygonMesh& pmesh, + const NamedParameters& np) { using boost::choose_param; using boost::get_param; @@ -153,7 +146,7 @@ compute_face_normal(typename boost::graph_traits::face_descriptor f Vector normal = traits.construct_vector_3_object()(CGAL::NULL_VECTOR); sum_normals(pmesh, f, vpmap, normal, traits); - if (!traits.equal_3_object()(normal, CGAL::NULL_VECTOR)) + if(!traits.equal_3_object()(normal, CGAL::NULL_VECTOR)) internal::normalize(normal, traits); return normal; @@ -167,8 +160,7 @@ typename GetGeomTraits::type::Vector_3 compute_face_normal(typename boost::graph_traits::face_descriptor f, const PolygonMesh& pmesh) { - return compute_face_normal(f, pmesh, - CGAL::Polygon_mesh_processing::parameters::all_default()); + return compute_face_normal(f, pmesh, CGAL::Polygon_mesh_processing::parameters::all_default()); } /// \endcond @@ -196,17 +188,15 @@ compute_face_normal(typename boost::graph_traits::face_descriptor f * If `Kernel::FT` does not have a `sqrt()` operation, the square root computation * will be done approximately. */ -template -void -compute_face_normals(const PolygonMesh& pmesh - , FaceNormalMap fnm - , const NamedParameters& np) +template +void compute_face_normals(const PolygonMesh& pmesh, + FaceNormalMap fnm, + const NamedParameters& np) { typedef typename GetGeomTraits::type Kernel; - for(typename boost::graph_traits::face_descriptor f : faces(pmesh)){ + for(typename boost::graph_traits::face_descriptor f : faces(pmesh)) + { typename Kernel::Vector_3 vec = compute_face_normal(f, pmesh, np); put(fnm, f, vec); } @@ -226,17 +216,16 @@ void compute_face_normals(const PolygonMesh& pmesh, FaceNormalMap fnm) template bool almost_equal(const typename GT::Vector_3 v1, const typename GT::Vector_3 v2) { - return (CGAL::abs(1 - GT().compute_scalar_product_3_object()(v1, v2)) < 1e-10); + return (CGAL::abs(1 - GT().compute_scalar_product_3_object()(v1, v2)) < 1e-10); // @tolerance } template -bool -does_enclose_other_normals(const int i, const int j, const int k, - const typename K::Vector_3& nb, - const typename K::FT sp_bi, - std::vector::face_descriptor> incident_faces, - const FaceNormalVector& normalized_face_normals, - const K& traits) +bool does_enclose_other_normals(const int i, const int j, const int k, + const typename K::Vector_3& nb, + const typename K::FT sp_bi, + std::vector::face_descriptor> incident_faces, + const FaceNormalVector& normalized_face_normals, + const K& traits) { typedef typename K::FT FT; typedef typename K::Vector_3 Vector_3; @@ -490,7 +479,8 @@ void compute_most_visible_vertex_normals(VertexNormalMap& vertex_normal_map, } template -void compute_most_visible_vertex_normals(VertexNormalMap& vertex_normal_map, const PolygonMesh& pmesh) +void compute_most_visible_vertex_normals(VertexNormalMap& vertex_normal_map, + const PolygonMesh& pmesh) { return compute_most_visible_vertex_normals(vertex_normal_map, pmesh, CGAL::parameters::all_default()); } @@ -528,8 +518,7 @@ typename GetGeomTraits::type::Vector_3 #endif compute_vertex_normal(typename boost::graph_traits::vertex_descriptor v, const PolygonMesh& pmesh, - const NamedParameters& np - ) + const NamedParameters& np) { using boost::choose_param; @@ -543,10 +532,7 @@ compute_vertex_normal(typename boost::graph_traits::vertex_descript NamedParameters, DefaultMap> ::type FaceNormalMap; FaceNormalMap fnmap = choose_param(get_param(np, internal_np::face_normal), DefaultMap()); - bool fnmap_valid - = !boost::is_same::value; + bool fnmap_valid = !boost::is_same::value; typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; @@ -554,11 +540,11 @@ compute_vertex_normal(typename boost::graph_traits::vertex_descript halfedge_descriptor he = halfedge(v, pmesh); // handle isolated vertices - if (he==boost::graph_traits::null_halfedge()) return normal; + if(he==boost::graph_traits::null_halfedge()) return normal; halfedge_descriptor end = he; do { - if (!is_border(he, pmesh)) + if(!is_border(he, pmesh)) { Vector n = fnmap_valid ? get(fnmap, face(he, pmesh)) : compute_face_normal(face(he, pmesh), pmesh, np); @@ -567,8 +553,9 @@ compute_vertex_normal(typename boost::graph_traits::vertex_descript he = opposite(next(he, pmesh), pmesh); } while (he != end); - if ( ! traits.equal_3_object()(normal, CGAL::NULL_VECTOR)) + if(!traits.equal_3_object()(normal, CGAL::NULL_VECTOR)) internal::normalize(normal, traits); + return normal; } @@ -608,19 +595,15 @@ compute_vertex_normal(typename boost::graph_traits::vertex_descript * If `Kernel::FT` does not have a `sqrt()` operation, the square root computation * will be done approximately. */ -template -void -compute_vertex_normals(const PolygonMesh& pmesh - , VertexNormalMap vnm - , const NamedParameters& np - ) +template +void compute_vertex_normals(const PolygonMesh& pmesh, + VertexNormalMap vnm, + const NamedParameters& np) { typedef typename GetGeomTraits::type Kernel; - for(typename boost::graph_traits::vertex_descriptor v : vertices(pmesh)){ + for(typename boost::graph_traits::vertex_descriptor v : vertices(pmesh)) + { typename Kernel::Vector_3 vec = compute_vertex_normal(v, pmesh, np); put(vnm, v, vec); } @@ -664,17 +647,11 @@ void compute_vertex_normals(const PolygonMesh& pmesh, VertexNormalMap vnm) * If `Kernel::FT` does not have a `sqrt()` operation, the square root computation * will be done approximately. */ -template -void -compute_normals(const PolygonMesh& pmesh - , VertexNormalMap vnm - , FaceNormalMap fnm - , const NamedParameters& np - ) +template +void compute_normals(const PolygonMesh& pmesh, + VertexNormalMap vnm, + FaceNormalMap fnm, + const NamedParameters& np) { compute_face_normals(pmesh, fnm, np); compute_vertex_normals(pmesh, vnm, np.face_normal_map(fnm)); From 3fe52056cfd5b025768156df43ae6b6570319534 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Thu, 18 Jul 2019 11:19:03 +0200 Subject: [PATCH 06/39] Do not normalize bisector in most visible vertex computations --- .../include/CGAL/Polygon_mesh_processing/compute_normal.h | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h index 59e02ff731f..782c8251e7c 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h @@ -253,15 +253,14 @@ bool does_enclose_other_normals(const int i, const int j, const int k, } template -typename GT::Vector_3 compute_normals_bisector(const typename GT::Vector_3 ni, - const typename GT::Vector_3 nj, +typename GT::Vector_3 compute_normals_bisector(const typename GT::Vector_3& ni, + const typename GT::Vector_3& nj, const GT& traits) { if(ni == nj) return ni; typename GT::Vector_3 nb = traits.construct_sum_of_vectors_3_object()(ni, nj); - CGAL::Polygon_mesh_processing::internal::normalize(nb, traits); // prob not necessary @fixme return nb; } From 52ca4f828f201778b321f65fbb1da01faf2b7985 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Mon, 12 Aug 2019 07:49:12 -0700 Subject: [PATCH 07/39] Comment debug output as travis fails --- .../include/CGAL/Polygon_mesh_processing/compute_normal.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h index 782c8251e7c..b68e3988e25 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h @@ -308,7 +308,7 @@ compute_most_visible_normal_2_points(std::vector typename GT::Vector_3 -compute_most_visible_normal_3_points(std::vector::face_descriptor> incident_faces, +compute_most_visible_normal_3_points(const std::vector::face_descriptor>& incident_faces, const FaceNormalVector& normalized_face_normals, const GT& traits) { @@ -463,7 +463,7 @@ void compute_most_visible_vertex_normals(VertexNormalMap& vertex_normal_map, normalized_face_normals[fd] = en; } - std::ofstream out("computed_normals.cgal.polylines.txt"); + // std::ofstream out("computed_normals.cgal.polylines.txt"); // we start by evaluating the translation normal for each vertex for(vertex_descriptor vd : vertices(pmesh)) @@ -473,7 +473,7 @@ void compute_most_visible_vertex_normals(VertexNormalMap& vertex_normal_map, CGAL::Polygon_mesh_processing::internal::normalize(n, traits); vertex_normal_map[vd] = n; - out << "2 " << pmesh.point(vd) << " " << pmesh.point(vd) + n << std::endl; + // out << "2 " << pmesh.point(vd) << " " << pmesh.point(vd) + n << std::endl; } } From 7db6c8e67088a4015dbd55304673a80cfbc3e39f Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Mon, 12 Aug 2019 07:51:32 -0700 Subject: [PATCH 08/39] make const & --- .../include/CGAL/Polygon_mesh_processing/compute_normal.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h index b68e3988e25..99e662ed4f5 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h @@ -214,7 +214,7 @@ void compute_face_normals(const PolygonMesh& pmesh, FaceNormalMap fnm) /// \endcond template -bool almost_equal(const typename GT::Vector_3 v1, const typename GT::Vector_3 v2) +bool almost_equal(const typename GT::Vector_3& v1, const typename GT::Vector_3& v2) { return (CGAL::abs(1 - GT().compute_scalar_product_3_object()(v1, v2)) < 1e-10); // @tolerance } @@ -223,7 +223,7 @@ template bool does_enclose_other_normals(const int i, const int j, const int k, const typename K::Vector_3& nb, const typename K::FT sp_bi, - std::vector::face_descriptor> incident_faces, + const std::vector::face_descriptor>& incident_faces, const FaceNormalVector& normalized_face_normals, const K& traits) { From f05a11dde52ed54659a8b689d06e8b04e17fc9a2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Fri, 23 Aug 2019 11:15:41 +0200 Subject: [PATCH 09/39] Fix some internal functions not being in the 'internal' namespace --- .../CGAL/Polygon_mesh_processing/compute_normal.h | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h index 747c0a6b184..27633df4f0f 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h @@ -66,8 +66,6 @@ triangle_normal(const Point& p0, const Point& p1, const Point& p2, const GT& tra return traits.construct_scaled_vector_3_object()(n, 0.5); } -} // namespace internal - template void sum_normals(const PM& pmesh, typename boost::graph_traits::face_descriptor f, @@ -93,6 +91,8 @@ void sum_normals(const PM& pmesh, } } +} // namespace internal + /** * \ingroup PMP_normal_grp * computes the outward unit vector normal to face `f`. @@ -143,7 +143,7 @@ compute_face_normal(typename boost::graph_traits::face_descriptor f typedef typename GT::Vector_3 Vector; Vector normal = traits.construct_vector_3_object()(CGAL::NULL_VECTOR); - sum_normals(pmesh, f, vpmap, normal, traits); + internal::sum_normals(pmesh, f, vpmap, normal, traits); if(!traits.equal_3_object()(normal, CGAL::NULL_VECTOR)) internal::normalize(normal, traits); @@ -202,6 +202,8 @@ void compute_face_normals(const PolygonMesh& pmesh, FaceNormalMap fnm) compute_face_normals(pmesh, fnm, CGAL::Polygon_mesh_processing::parameters::all_default()); } +namespace internal { + template bool almost_equal(const typename GT::Vector_3& v1, const typename GT::Vector_3& v2) { @@ -466,12 +468,7 @@ void compute_most_visible_vertex_normals(VertexNormalMap& vertex_normal_map, } } -template -void compute_most_visible_vertex_normals(VertexNormalMap& vertex_normal_map, - const PolygonMesh& pmesh) -{ - return compute_most_visible_vertex_normals(vertex_normal_map, pmesh, CGAL::parameters::all_default()); -} +} // end namespace internal /** * \ingroup PMP_normal_grp From ad90fa30fd064876404dcce283684e711540b0ba Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Fri, 23 Aug 2019 11:16:01 +0200 Subject: [PATCH 10/39] Fix potentially taking references to temporaries For example, if the vertex point property map returned temporaries --- .../include/CGAL/Polygon_mesh_processing/compute_normal.h | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h index 27633df4f0f..31b8b158a24 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h @@ -76,13 +76,15 @@ void sum_normals(const PM& pmesh, typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + typedef boost::property_traits::reference Point_ref; + halfedge_descriptor he = halfedge(f, pmesh); vertex_descriptor v = source(he, pmesh); - const Point& pv = get(vpmap, v); + const Point_ref pv = get(vpmap, v); while(v != target(next(he, pmesh), pmesh)) { - const Point& pvn = get(vpmap, target(he, pmesh)); - const Point& pvnn = get(vpmap, target(next(he, pmesh), pmesh)); + const Point_ref pvn = get(vpmap, target(he, pmesh)); + const Point_ref pvnn = get(vpmap, target(next(he, pmesh), pmesh)); Vector n = internal::triangle_normal(pv, pvn, pvnn, traits); sum = traits.construct_sum_of_vectors_3_object()(sum, n); From bba6e57053e85fef7604638fa142990f893aedee Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Fri, 23 Aug 2019 11:35:34 +0200 Subject: [PATCH 11/39] Misc code cleaning --- .../Polygon_mesh_processing/compute_normal.h | 101 ++++++++++-------- 1 file changed, 55 insertions(+), 46 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h index 31b8b158a24..c3eb7ec4a50 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h @@ -1,4 +1,4 @@ -// Copyright (c) 2015 GeometryFactory (France). +// Copyright (c) 2015-2019 GeometryFactory (France). // All rights reserved. // // This file is part of CGAL (www.cgal.org). @@ -15,26 +15,28 @@ // $URL$ // $Id$ // SPDX-License-Identifier: GPL-3.0+ -// // // Author(s) : Jane Tournois +// Mael Rouxel-Labbé #ifndef CGAL_POLYGON_MESH_PROCESSING_COMPUTE_NORMAL_H #define CGAL_POLYGON_MESH_PROCESSING_COMPUTE_NORMAL_H #include -#include -#include -#include -#include -#include -#include - #include #include -#include +#include +#include +#include + +#include + +#include +#include +#include +#include namespace CGAL { namespace Polygon_mesh_processing { @@ -43,10 +45,12 @@ namespace internal { template void normalize(typename GT::Vector_3& v, const GT& traits) { - typename GT::FT norm = CGAL::approximate_sqrt(traits.compute_squared_length_3_object()(v)); + typedef typename GT::FT FT; + //If the vector is small enough, approx_sqrt might return 0, and then we get nan values. //To avoid that, we check the resulted norm. If it is 0, we don't normalize. - if(norm != 0) + const FT norm = CGAL::approximate_sqrt(traits.compute_squared_length_3_object()(v)); + if(norm != FT(0)) { v = traits.construct_divided_vector_3_object()(v, norm); } @@ -56,6 +60,8 @@ template typename GT::Vector_3 triangle_normal(const Point& p0, const Point& p1, const Point& p2, const GT& traits) { + typedef typename GT::FT FT; + typename GT::Vector_3 n = traits.construct_cross_product_vector_3_object()( traits.construct_vector_3_object()(p1, p2), traits.construct_vector_3_object()(p1, p0)); @@ -63,7 +69,7 @@ triangle_normal(const Point& p0, const Point& p1, const Point& p2, const GT& tra //cross-product(AB, AC)'s norm is the area of the parallelogram //formed by these 2 vectors. //the triangle's area is half of it - return traits.construct_scaled_vector_3_object()(n, 0.5); + return traits.construct_scaled_vector_3_object()(n, FT(1)/FT(2)); } template @@ -207,9 +213,10 @@ void compute_face_normals(const PolygonMesh& pmesh, FaceNormalMap fnm) namespace internal { template -bool almost_equal(const typename GT::Vector_3& v1, const typename GT::Vector_3& v2) +bool almost_equal(const typename GT::Vector_3& v1, const typename GT::Vector_3& v2, + const GT& traits) { - return (CGAL::abs(1 - GT().compute_scalar_product_3_object()(v1, v2)) < 1e-10); // @tolerance + return (CGAL::abs(1 - traits.compute_scalar_product_3_object()(v1, v2)) < 1e-10); // @tolerance } template @@ -330,18 +337,18 @@ compute_most_visible_normal_3_points(const std::vector(ni, nj)) + if(almost_equal(ni, nj, traits)) { - if(almost_equal(nj, nk)) + if(almost_equal(nj, nk, traits)) nb = ni; else // ni == nj, but nj != nk nb = compute_normals_bisector(nj, nk, traits); } - else if(almost_equal(ni, nk)) // ni != nj + else if(almost_equal(ni, nk, traits)) // ni != nj { nb = compute_normals_bisector(nj, nk, traits); } - else if(almost_equal(nj, nk)) // ni != nj + else if(almost_equal(nj, nk, traits)) // ni != nj { nb = compute_normals_bisector(ni, nk, traits); } @@ -350,9 +357,13 @@ compute_most_visible_normal_3_points(const std::vector +template typename GT::Vector_3 compute_vertex_normal_most_visible_min_circle(typename boost::graph_traits::vertex_descriptor vd, const FaceNormalVector& normalized_face_normals, @@ -404,7 +415,7 @@ compute_vertex_normal_most_visible_min_circle(typename boost::graph_traits incident_faces; for(face_descriptor fd : CGAL::faces_around_target(hd, pmesh)) @@ -437,21 +448,17 @@ void compute_most_visible_vertex_normals(VertexNormalMap& vertex_normal_map, CGAL_precondition(CGAL::is_triangle_mesh(pmesh)); - typedef typename GetVertexPointMap::const_type VPMap; -// VPMap vpmap = choose_param(get_param(np, internal_np::vertex_point), -// get_const_property_map(vertex_point, pmesh)); - typedef typename GetGeomTraits::type GT; const GT traits = choose_param(get_param(np, internal_np::geom_traits), GT()); - typedef typename GT::Vector_3 Vector; + typedef typename GT::Vector_3 Vector_3; std::cout << "compute face normals" << std::endl; - std::vector normalized_face_normals(num_faces(pmesh)); + std::vector normalized_face_normals(num_faces(pmesh)); for(face_descriptor fd : faces(pmesh)) { - Vector en = compute_face_normal(fd, pmesh); + Vector en = compute_face_normal(fd, pmesh, np); CGAL::Polygon_mesh_processing::internal::normalize(en, traits); normalized_face_normals[fd] = en; } @@ -461,8 +468,7 @@ void compute_most_visible_vertex_normals(VertexNormalMap& vertex_normal_map, // we start by evaluating the translation normal for each vertex for(vertex_descriptor vd : vertices(pmesh)) { - Vector n = compute_vertex_normal_most_visible_min_circle(vd, normalized_face_normals, pmesh, traits); - + Vector_3 n = compute_vertex_normal_most_visible_min_circle(vd, normalized_face_normals, pmesh, traits); CGAL::Polygon_mesh_processing::internal::normalize(n, traits); vertex_normal_map[vd] = n; @@ -510,25 +516,26 @@ compute_vertex_normal(typename boost::graph_traits::vertex_descript using parameters::choose_parameter; using parameters::get_parameter; - typedef typename GetGeomTraits::type GT; - typedef typename GT::Vector_3 Vector; + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + + typedef typename GetGeomTraits::type GT; + typedef typename GT::Vector_3 Vector; GT traits = choose_parameter(get_parameter(np, internal_np::geom_traits), GT()); - typedef typename GetFaceNormalMap::NoMap DefaultMap; - typedef typename internal_np::Lookup_named_param_def < - internal_np::face_normal_t, - NamedParameters, - DefaultMap> ::type FaceNormalMap; + typedef typename GetFaceNormalMap::NoMap DefaultMap; + typedef typename internal_np::Lookup_named_param_def ::type FaceNormalMap; FaceNormalMap fnmap = choose_parameter(get_parameter(np, internal_np::face_normal), DefaultMap()); - bool fnmap_valid = !boost::is_same::value; - - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + const bool fnmap_valid = !std::is_same::value; Vector normal = traits.construct_vector_3_object()(CGAL::NULL_VECTOR); halfedge_descriptor he = halfedge(v, pmesh); // handle isolated vertices - if(he==boost::graph_traits::null_halfedge()) return normal; + if(he == boost::graph_traits::null_halfedge()) + return normal; + halfedge_descriptor end = he; do { @@ -538,8 +545,10 @@ compute_vertex_normal(typename boost::graph_traits::vertex_descript : compute_face_normal(face(he, pmesh), pmesh, np); normal = traits.construct_sum_of_vectors_3_object()(normal, n); } + he = opposite(next(he, pmesh), pmesh); - } while (he != end); + } + while (he != end); if(!traits.equal_3_object()(normal, CGAL::NULL_VECTOR)) internal::normalize(normal, traits); From cf357713371ca61bebc2ea650d5e65c4ca66474d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Fri, 23 Aug 2019 16:10:38 +0200 Subject: [PATCH 12/39] Robustify and incorporate most visible normal computations into standard fns --- .../Polygon_mesh_processing/compute_normal.h | 416 +++++++++++------- 1 file changed, 253 insertions(+), 163 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h index c3eb7ec4a50..4d52a94160f 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h @@ -29,15 +29,19 @@ #include #include +#include #include #include #include -#include #include #include +// @tmp +#include +#include + namespace CGAL { namespace Polygon_mesh_processing { namespace internal { @@ -79,10 +83,10 @@ void sum_normals(const PM& pmesh, Vector& sum, const GT& traits) { - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - typedef boost::property_traits::reference Point_ref; + typedef typename boost::property_traits::reference Point_ref; halfedge_descriptor he = halfedge(f, pmesh); vertex_descriptor v = source(he, pmesh); @@ -92,7 +96,7 @@ void sum_normals(const PM& pmesh, const Point_ref pvn = get(vpmap, target(he, pmesh)); const Point_ref pvnn = get(vpmap, target(next(he, pmesh), pmesh)); - Vector n = internal::triangle_normal(pv, pvn, pvnn, traits); + const Vector n = internal::triangle_normal(pv, pvn, pvnn, traits); sum = traits.construct_sum_of_vectors_3_object()(sum, n); he = next(he, pmesh); @@ -140,17 +144,17 @@ compute_face_normal(typename boost::graph_traits::face_descriptor f using parameters::choose_parameter; using parameters::get_parameter; - typedef typename GetGeomTraits::type GT; + typedef typename GetGeomTraits::type GT; GT traits = choose_parameter(get_parameter(np, internal_np::geom_traits), GT()); - typedef typename GetVertexPointMap::const_type VPMap; + typedef typename GetVertexPointMap::const_type VPMap; VPMap vpmap = choose_parameter(get_parameter(np, internal_np::vertex_point), get_const_property_map(vertex_point, pmesh)); - typedef typename GT::Point_3 Point; - typedef typename GT::Vector_3 Vector; + typedef typename GT::Point_3 Point; + typedef typename GT::Vector_3 Vector_3; - Vector normal = traits.construct_vector_3_object()(CGAL::NULL_VECTOR); + Vector_3 normal = traits.construct_vector_3_object()(CGAL::NULL_VECTOR); internal::sum_normals(pmesh, f, vpmap, normal, traits); if(!traits.equal_3_object()(normal, CGAL::NULL_VECTOR)) @@ -171,12 +175,12 @@ compute_face_normal(typename boost::graph_traits::face_descriptor f * \ingroup PMP_normal_grp * computes the outward unit vector normal for all faces of the polygon mesh. * @tparam PolygonMesh a model of `FaceGraph` -* @tparam FaceNormalMap a model of `WritablePropertyMap` with +* @tparam Face_normal_map a model of `WritablePropertyMap` with `boost::graph_traits::%face_descriptor` as key type and `Kernel::Vector_3` as value type. * * @param pmesh the polygon mesh -* @param fnm the property map in which the normals are written +* @param face_normals the property map in which the normals are written * @param np optional sequence of \ref pmp_namedparameters "Named Parameters" among the ones listed below * * \cgalNamedParamsBegin @@ -190,9 +194,9 @@ compute_face_normal(typename boost::graph_traits::face_descriptor f * If `Kernel::FT` does not have a `sqrt()` operation, the square root computation * will be done approximately. */ -template +template void compute_face_normals(const PolygonMesh& pmesh, - FaceNormalMap fnm, + Face_normal_map face_normals, const NamedParameters& np) { typedef typename GetGeomTraits::type Kernel; @@ -200,14 +204,14 @@ void compute_face_normals(const PolygonMesh& pmesh, for(typename boost::graph_traits::face_descriptor f : faces(pmesh)) { typename Kernel::Vector_3 vec = compute_face_normal(f, pmesh, np); - put(fnm, f, vec); + put(face_normals, f, vec); } } -template -void compute_face_normals(const PolygonMesh& pmesh, FaceNormalMap fnm) +template +void compute_face_normals(const PolygonMesh& pmesh, Face_normal_map face_normals) { - compute_face_normals(pmesh, fnm, CGAL::Polygon_mesh_processing::parameters::all_default()); + compute_face_normals(pmesh, face_normals, CGAL::Polygon_mesh_processing::parameters::all_default()); } namespace internal { @@ -216,7 +220,16 @@ template bool almost_equal(const typename GT::Vector_3& v1, const typename GT::Vector_3& v2, const GT& traits) { - return (CGAL::abs(1 - traits.compute_scalar_product_3_object()(v1, v2)) < 1e-10); // @tolerance + typedef typename GT::FT FT; + + // We are doing a lot of likely inexact constructions to compute min circles on the sphere and + // these degenerate cases often happen (e.g. almost flat surfaces), so we don't want to to the usual + // "switch to the exact kernel" in degenerate cases robustness trick. Rather, use a fixed tolerance + // to assimilate vectors and scalar products. + // + // Tolerance is (arbitrarily) chosen as theta := 0.01°, thus we want sp(v1,v2) >= cos(theta) + const FT cos_theta = 0.99999998476912910; + return traits.compute_scalar_product_3_object()(v1, v2) >= cos_theta; } template @@ -224,25 +237,33 @@ bool does_enclose_other_normals(const int i, const int j, const int k, const typename K::Vector_3& nb, const typename K::FT sp_bi, const std::vector::face_descriptor>& incident_faces, - const FaceNormalVector& normalized_face_normals, + const FaceNormalVector& face_normals, const K& traits) { - typedef typename K::FT FT; - typedef typename K::Vector_3 Vector_3; + typedef typename K::FT FT; + typedef typename K::Vector_3 Vector_3; + typedef typename boost::property_traits::reference Vector_ref; typename K::Compute_scalar_product_3 sp = traits.compute_scalar_product_3_object(); + const FT nbn = CGAL::sqrt(traits.compute_squared_length_3_object()(nb)); + // check that this min circle defined by the diameter contains the other points const std::size_t nif = incident_faces.size(); for(int l=0; l +typename GT::Vector_3 compute_normals_bisector(const typename GT::Vector_3& ni, + const typename GT::Vector_3& nj, + const typename GT::Vector_3& nk, + const GT& traits) +{ + typedef typename GT::Point_3 Point_3; + typedef typename GT::Vector_3 Vector_3; + + Vector_3 nb = traits.construct_vector_3_object()(CGAL::NULL_VECTOR); + + if(almost_equal(ni, nj, traits)) + { + if(almost_equal(nj, nk, traits)) + nb = ni; + else // ni == nj, but nij != nk + nb = compute_normals_bisector(nj, nk, traits); + } + else if(almost_equal(ni, nk, traits)) // ni != nj + { + nb = compute_normals_bisector(nj, nk, traits); + } + else if(almost_equal(nj, nk, traits)) // ni != nj, ni != nk + { + nb = compute_normals_bisector(ni, nk, traits); + } + else + { + CGAL_assertion(ni != nj); + CGAL_assertion(ni != nk); + CGAL_assertion(nj != nk); + CGAL_assertion(!traits.collinear_3_object()(CGAL::ORIGIN + ni, + CGAL::ORIGIN + nj, + CGAL::ORIGIN + nk)); + +#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP + std::cout << "Triplet: ni[" << ni << "] nj[" << nj << "] nk[" << nk << "]" << std::endl; +#endif + + const Point_3 c = traits.construct_circumcenter_3_object()(CGAL::ORIGIN + ni, + CGAL::ORIGIN + nj, + CGAL::ORIGIN + nk); + CGAL_assertion(c != CGAL::ORIGIN); + + nb = traits.construct_vector_3_object()(CGAL::ORIGIN, c); // note that this isn't normalized + } return nb; } template std::pair -compute_most_visible_normal_2_points(std::vector::face_descriptor> incident_faces, - const FaceNormalVector& normalized_face_normals, +compute_most_visible_normal_2_points(std::vector::face_descriptor>& incident_faces, + const FaceNormalVector& face_normals, const GT& traits) { typedef typename GT::FT FT; typedef typename GT::Vector_3 Vector_3; + typedef typename boost::property_traits::reference Vector_ref; typename GT::Compute_scalar_product_3 sp = traits.compute_scalar_product_3_object(); - FT min_sp = -1; +#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP + std::cout << "Trying to find enclosing normal with 2 normals" << std::endl; +#endif - Vector_3 n = CGAL::NULL_VECTOR; + FT min_sp = -1; + Vector_3 n = traits.construct_vector_3_object()(CGAL::NULL_VECTOR); const std::size_t nif = incident_faces.size(); for(int i=0; i= 0); if(sp_bi <= min_sp) continue; - if(!does_enclose_other_normals(i, j, -1 /*NA*/, nb, sp_bi, incident_faces, normalized_face_normals, traits)) + if(!does_enclose_other_normals(i, j, -1 /*NA*/, nb, sp_bi, incident_faces, face_normals, traits)) continue; min_sp = sp_bi; @@ -303,25 +383,30 @@ compute_most_visible_normal_2_points(std::vector typename GT::Vector_3 compute_most_visible_normal_3_points(const std::vector::face_descriptor>& incident_faces, - const FaceNormalVector& normalized_face_normals, + const FaceNormalVector& face_normals, const GT& traits) { typedef typename GT::FT FT; - typedef typename GT::Point_3 Point_3; typedef typename GT::Vector_3 Vector_3; + typedef typename boost::property_traits::reference Vector_ref; typename GT::Compute_scalar_product_3 sp = traits.compute_scalar_product_3_object(); typename GT::Construct_sum_of_vectors_3 csv = traits.construct_sum_of_vectors_3_object(); +#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP + std::cout << "Trying to find enclosing normal with 3 normals" << std::endl; +#endif + FT min_sp = -1; - Vector_3 n = CGAL::NULL_VECTOR; + Vector_3 n = traits.construct_vector_3_object()(CGAL::NULL_VECTOR); const std::size_t nif = incident_faces.size(); @@ -331,48 +416,15 @@ compute_most_visible_normal_3_points(const std::vector(i, j, k, nb, sp_bi, incident_faces, normalized_face_normals, traits)) + if(!does_enclose_other_normals(i, j, k, nb, sp_bi, incident_faces, face_normals, traits)) continue; min_sp = sp_bi; @@ -391,89 +443,73 @@ compute_most_visible_normal_3_points(const std::vector typename GT::Vector_3 -compute_vertex_normal_most_visible_min_circle(typename boost::graph_traits::vertex_descriptor vd, - const FaceNormalVector& normalized_face_normals, +compute_vertex_normal_most_visible_min_circle(typename boost::graph_traits::vertex_descriptor v, + const FaceNormalVector& face_normals, const PolygonMesh& pmesh, const GT& traits) { + typedef typename boost::graph_traits::face_descriptor face_descriptor; + typedef typename GT::FT FT; typedef typename GT::Vector_3 Vector_3; - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - typedef typename boost::graph_traits::face_descriptor face_descriptor; - -#ifdef CGAL_MOST_VISIBLE_NORMAL_VERBOSE - std::cout << std::endl << std::endl; - std::cout << "----------------------------------------------------------------------" << std::endl; - std::cout << "compute_vertex_normal_most_visible_min_circle at " << pmesh.point(vd) << std::endl; -#endif - - const halfedge_descriptor hd = halfedge(vd, pmesh); - std::vector incident_faces; - for(face_descriptor fd : CGAL::faces_around_target(hd, pmesh)) + for(face_descriptor f : CGAL::faces_around_target(halfedge(v, pmesh), pmesh)) { - if(fd == boost::graph_traits::null_face()) + if(f == boost::graph_traits::null_face()) continue; - incident_faces.push_back(fd); + incident_faces.push_back(f); } if(incident_faces.size() == 1) - return normalized_face_normals[incident_faces.front()]; + return get(face_normals, incident_faces.front()); - std::pair res = compute_most_visible_normal_2_points(incident_faces, normalized_face_normals, traits); - if(res.second) + std::pair res = compute_most_visible_normal_2_points(incident_faces, face_normals, traits); + + if(!res.second || // won't be able to find a most visible normal + (res.second && res.first != CGAL::NULL_VECTOR)) // found a valid normal through 2 point min circle return res.first; CGAL_assertion(incident_faces.size() > 2); - return compute_most_visible_normal_3_points(incident_faces, normalized_face_normals, traits); + return compute_most_visible_normal_3_points(incident_faces, face_normals, traits); } -template -void compute_most_visible_vertex_normals(VertexNormalMap& vertex_normal_map, - const PolygonMesh& pmesh, - const NamedParameters& np) +template +typename GT::Vector_3 +compute_vertex_normal_as_weighted_sum_of_face_normals(typename boost::graph_traits::vertex_descriptor v, + const FaceNormalVector& face_normals, + const PolygonMesh& pmesh, + const GT& traits) { - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - typedef typename boost::graph_traits::face_descriptor face_descriptor; + typedef typename boost::graph_traits::face_descriptor face_descriptor; - CGAL_precondition(CGAL::is_triangle_mesh(pmesh)); + typedef typename GT::Vector_3 Vector_3; + typedef typename boost::property_traits::reference Vector_ref; - typedef typename GetGeomTraits::type GT; - const GT traits = choose_param(get_param(np, internal_np::geom_traits), GT()); + Vector_3 normal = traits.construct_vector_3_object()(CGAL::NULL_VECTOR); - typedef typename GT::Vector_3 Vector_3; - - std::cout << "compute face normals" << std::endl; - - std::vector normalized_face_normals(num_faces(pmesh)); - for(face_descriptor fd : faces(pmesh)) + for(face_descriptor f : CGAL::faces_around_target(halfedge(v, pmesh), pmesh)) { - Vector en = compute_face_normal(fd, pmesh, np); - CGAL::Polygon_mesh_processing::internal::normalize(en, traits); - normalized_face_normals[fd] = en; + if(f == boost::graph_traits::null_face()) + continue; + + const Vector_ref n = get(face_normals, f); + normal = traits.construct_sum_of_vectors_3_object()(normal, n); } - // std::ofstream out("computed_normals.cgal.polylines.txt"); - - // we start by evaluating the translation normal for each vertex - for(vertex_descriptor vd : vertices(pmesh)) - { - Vector_3 n = compute_vertex_normal_most_visible_min_circle(vd, normalized_face_normals, pmesh, traits); - CGAL::Polygon_mesh_processing::internal::normalize(n, traits); - vertex_normal_map[vd] = n; - - // out << "2 " << pmesh.point(vd) << " " << pmesh.point(vd) + n << std::endl; - } + return normal; } } // end namespace internal @@ -514,41 +550,64 @@ compute_vertex_normal(typename boost::graph_traits::vertex_descript const NamedParameters& np) { using parameters::choose_parameter; + using parameters::is_default_parameter; using parameters::get_parameter; typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + typedef typename boost::graph_traits::face_descriptor face_descriptor; typedef typename GetGeomTraits::type GT; - typedef typename GT::Vector_3 Vector; + typedef typename GT::Vector_3 Vector_3; GT traits = choose_parameter(get_parameter(np, internal_np::geom_traits), GT()); - typedef typename GetFaceNormalMap::NoMap DefaultMap; + typedef typename GetVertexPointMap::const_type VPMap; + VPMap vpmap = choose_parameter(get_parameter(np, internal_np::vertex_point), + get_const_property_map(vertex_point, pmesh)); + + typedef std::map Face_vector_map; + typedef boost::associative_property_map Default_map; + typedef typename internal_np::Lookup_named_param_def ::type FaceNormalMap; - FaceNormalMap fnmap = choose_parameter(get_parameter(np, internal_np::face_normal), DefaultMap()); - const bool fnmap_valid = !std::is_same::value; + Default_map>::type Face_normal_map; + Face_normal_map face_normals = choose_parameter(get_parameter(np, internal_np::face_normal), + Default_map()); + const bool fnmap_valid = is_default_parameter(get_parameter(np, internal_np::face_normal)); - Vector normal = traits.construct_vector_3_object()(CGAL::NULL_VECTOR); +#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP + std::cout << std::endl << std::endl; + std::cout << "----------------------------------------------------------------------" << std::endl; + std::cout << "compute vertex at " << get(vpmap, v) << std::endl; +#endif - halfedge_descriptor he = halfedge(v, pmesh); // handle isolated vertices + halfedge_descriptor he = halfedge(v, pmesh); if(he == boost::graph_traits::null_halfedge()) - return normal; + return CGAL::NULL_VECTOR; - halfedge_descriptor end = he; - do + if(!fnmap_valid) { - if(!is_border(he, pmesh)) + for(face_descriptor f : CGAL::faces_around_target(halfedge(v, pmesh), pmesh)) { - Vector n = fnmap_valid ? get(fnmap, face(he, pmesh)) - : compute_face_normal(face(he, pmesh), pmesh, np); - normal = traits.construct_sum_of_vectors_3_object()(normal, n); - } + if(f == boost::graph_traits::null_face()) + continue; - he = opposite(next(he, pmesh), pmesh); + put(face_normals, f, compute_face_normal(f, pmesh, np)); + } } - while (he != end); + +#if 0 // old approach @tmp + Vector_3 normal = internal::compute_vertex_normal_as_weighted_sum_of_face_normals(v, face_normals, pmesh, traits); +#else + Vector_3 normal = internal::compute_vertex_normal_most_visible_min_circle(v, face_normals, pmesh, traits); + if(traits.equal_3_object()(normal, CGAL::NULL_VECTOR)) // can't always find a most visible normal + { +#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP + std::cout << "Failed to find most visible normal, use weighted sum of normals" << std::endl; +#endif + normal = internal::compute_vertex_normal_as_weighted_sum_of_face_normals(v, face_normals, pmesh, traits); + } +#endif if(!traits.equal_3_object()(normal, CGAL::NULL_VECTOR)) internal::normalize(normal, traits); @@ -589,22 +648,51 @@ compute_vertex_normal(typename boost::graph_traits::vertex_descript */ template void compute_vertex_normals(const PolygonMesh& pmesh, - VertexNormalMap vnm, + VertexNormalMap vertex_normals, const NamedParameters& np) { - typedef typename GetGeomTraits::type Kernel; + using parameters::choose_parameter; + using parameters::is_default_parameter; + using parameters::get_parameter; - for(typename boost::graph_traits::vertex_descriptor v : vertices(pmesh)) + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + + typedef typename GetGeomTraits::type GT; + typedef typename GT::Vector_3 Vector_3; + + typedef CGAL::dynamic_face_property_t Face_normal_tag; + typedef typename boost::property_map::type Face_normal_dmap; + + typedef typename internal_np::Lookup_named_param_def::type Face_normal_map; + Face_normal_map face_normals = choose_parameter(get_parameter(np, internal_np::face_normal), + get(Face_normal_tag(), pmesh)); + const bool fnmap_valid = is_default_parameter(get_parameter(np, internal_np::face_normal)); + + if(!fnmap_valid) + compute_face_normals(pmesh, face_normals, np); + + // @tmp ---- + std::ofstream out("computed_normals.cgal.polylines.txt"); + const Bbox_3 bb = bbox(pmesh, np); + const typename GT::FT bbox_diagonal = CGAL::sqrt(CGAL::square(bb.xmax() - bb.xmin()) + + CGAL::square(bb.ymax() - bb.ymin()) + + CGAL::square(bb.zmax() - bb.zmin())); + // ------ + + for(vertex_descriptor v : vertices(pmesh)) { - typename Kernel::Vector_3 vec = compute_vertex_normal(v, pmesh, np); - put(vnm, v, vec); + const Vector_3 n = compute_vertex_normal(v, pmesh, np.face_normal_map(face_normals)); + out << "2 " << pmesh.point(v) << " " << pmesh.point(v) + 0.1 * bbox_diagonal * n << std::endl; + put(vertex_normals, v, n); } } template -void compute_vertex_normals(const PolygonMesh& pmesh, VertexNormalMap vnm) +void compute_vertex_normals(const PolygonMesh& pmesh, VertexNormalMap vertex_normals) { - compute_vertex_normals(pmesh, vnm, CGAL::Polygon_mesh_processing::parameters::all_default()); + compute_vertex_normals(pmesh, vertex_normals, CGAL::parameters::all_default()); } /** @@ -619,8 +707,8 @@ void compute_vertex_normals(const PolygonMesh& pmesh, VertexNormalMap vnm) `Kernel::Vector_3` as value type. * * @param pmesh the polygon mesh -* @param vnm the property map in which the vertex normals are written -* @param fnm the property map in which the face normals are written +* @param vertex_normals the property map in which the vertex normals are written +* @param face_normals the property map in which the face normals are written * @param np optional sequence of \ref pmp_namedparameters "Named Parameters" among the ones listed below * * \cgalNamedParamsBegin @@ -634,22 +722,24 @@ void compute_vertex_normals(const PolygonMesh& pmesh, VertexNormalMap vnm) * If `Kernel::FT` does not have a `sqrt()` operation, the square root computation * will be done approximately. */ -template +template void compute_normals(const PolygonMesh& pmesh, - VertexNormalMap vnm, - FaceNormalMap fnm, + VertexNormalMap vertex_normals, + FaceNormalMap face_normals, const NamedParameters& np) { - compute_face_normals(pmesh, fnm, np); - compute_vertex_normals(pmesh, vnm, np.face_normal_map(fnm)); + compute_face_normals(pmesh, face_normals, np); + compute_vertex_normals(pmesh, vertex_normals, np.face_normal_map(face_normals)); } template void compute_normals(const PolygonMesh& pmesh, - VertexNormalMap vnm, - FaceNormalMap fnm) + VertexNormalMap vertex_normals, + FaceNormalMap face_normals) { - compute_normals(pmesh, vnm, fnm, CGAL::Polygon_mesh_processing::parameters::all_default()); + compute_normals(pmesh, vertex_normals, face_normals, CGAL::parameters::all_default()); } } // namespace Polygon_mesh_processing From 14e33d62ecfdcfb5b2579b2c5e6fe180859174c0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Mon, 26 Aug 2019 13:01:10 +0200 Subject: [PATCH 13/39] Add a convenience default template parameter --- BGL/include/CGAL/boost/graph/named_params_helper.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/BGL/include/CGAL/boost/graph/named_params_helper.h b/BGL/include/CGAL/boost/graph/named_params_helper.h index bf8c2c606da..ee8f5db1a0e 100644 --- a/BGL/include/CGAL/boost/graph/named_params_helper.h +++ b/BGL/include/CGAL/boost/graph/named_params_helper.h @@ -115,7 +115,8 @@ namespace CGAL { namespace Polygon_mesh_processing { - template + template > class GetVertexPointMap { typedef typename property_map_selector::const_type From b9f37cea38393796bb886ce31a9b9f613489bdeb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Mon, 26 Aug 2019 13:03:22 +0200 Subject: [PATCH 14/39] Hide debug code behind macros --- .../Polygon_mesh_processing/compute_normal.h | 23 ++++++++++++------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h index 4d52a94160f..92241aa1b08 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h @@ -38,9 +38,16 @@ #include #include -// @tmp +#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP +# ifndef CGAL_PMP_COMPUTE_NORMAL_DEBUG +# define CGAL_PMP_COMPUTE_NORMAL_DEBUG +# endif +#endif + +#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG #include #include +#endif namespace CGAL { namespace Polygon_mesh_processing { @@ -596,9 +603,6 @@ compute_vertex_normal(typename boost::graph_traits::vertex_descript } } -#if 0 // old approach @tmp - Vector_3 normal = internal::compute_vertex_normal_as_weighted_sum_of_face_normals(v, face_normals, pmesh, traits); -#else Vector_3 normal = internal::compute_vertex_normal_most_visible_min_circle(v, face_normals, pmesh, traits); if(traits.equal_3_object()(normal, CGAL::NULL_VECTOR)) // can't always find a most visible normal { @@ -607,7 +611,6 @@ compute_vertex_normal(typename boost::graph_traits::vertex_descript #endif normal = internal::compute_vertex_normal_as_weighted_sum_of_face_normals(v, face_normals, pmesh, traits); } -#endif if(!traits.equal_3_object()(normal, CGAL::NULL_VECTOR)) internal::normalize(normal, traits); @@ -673,19 +676,23 @@ void compute_vertex_normals(const PolygonMesh& pmesh, if(!fnmap_valid) compute_face_normals(pmesh, face_normals, np); - // @tmp ---- +#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG std::ofstream out("computed_normals.cgal.polylines.txt"); const Bbox_3 bb = bbox(pmesh, np); const typename GT::FT bbox_diagonal = CGAL::sqrt(CGAL::square(bb.xmax() - bb.xmin()) + CGAL::square(bb.ymax() - bb.ymin()) + CGAL::square(bb.zmax() - bb.zmin())); - // ------ +#endif for(vertex_descriptor v : vertices(pmesh)) { const Vector_3 n = compute_vertex_normal(v, pmesh, np.face_normal_map(face_normals)); - out << "2 " << pmesh.point(v) << " " << pmesh.point(v) + 0.1 * bbox_diagonal * n << std::endl; put(vertex_normals, v, n); + +#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG + out << "2 " << get(vpmap, v) << " " + << get(vpmap, v) + traits.construct_scaled_vector_3_object()(n, 0.1 * bbox_diagonal) << "\n"; +#endif } } From ce9b3a890d4f636aa44d71ff168c4afbb4559cbb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Mon, 26 Aug 2019 13:04:24 +0200 Subject: [PATCH 15/39] Use sin-based weights in the default vertex normal computation --- .../Polygon_mesh_processing/compute_normal.h | 68 +++++++++++++++---- 1 file changed, 54 insertions(+), 14 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h index 92241aa1b08..f73bb23492f 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h @@ -218,11 +218,17 @@ void compute_face_normals(const PolygonMesh& pmesh, template void compute_face_normals(const PolygonMesh& pmesh, Face_normal_map face_normals) { - compute_face_normals(pmesh, face_normals, CGAL::Polygon_mesh_processing::parameters::all_default()); + compute_face_normals(pmesh, face_normals, CGAL::parameters::all_default()); } namespace internal { +enum Vertex_normal_type { + NO_WEIGHT = 0, + SIN_WEIGHT, + MOST_VISIBLE +}; + template bool almost_equal(const typename GT::Vector_3& v1, const typename GT::Vector_3& v2, const GT& traits) @@ -493,28 +499,62 @@ compute_vertex_normal_most_visible_min_circle(typename boost::graph_traits(incident_faces, face_normals, traits); } -template +template typename GT::Vector_3 -compute_vertex_normal_as_weighted_sum_of_face_normals(typename boost::graph_traits::vertex_descriptor v, - const FaceNormalVector& face_normals, - const PolygonMesh& pmesh, - const GT& traits) +compute_vertex_normal_as_sum_of_weighted_normals(typename boost::graph_traits::vertex_descriptor v, + const Vertex_normal_type& vn_type, + const FaceNormalVector& face_normals, + const VertexPointMap& vpmap, + const PolygonMesh& pmesh, + const GT& traits) { - typedef typename boost::graph_traits::face_descriptor face_descriptor; + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + typedef typename GT::FT FT; typedef typename GT::Vector_3 Vector_3; typedef typename boost::property_traits::reference Vector_ref; - Vector_3 normal = traits.construct_vector_3_object()(CGAL::NULL_VECTOR); + typename GT::Construct_vector_3 cv_3 = traits.construct_vector_3_object(); + typename GT::Compute_squared_length_3 csl_3 = traits.compute_squared_length_3_object(); - for(face_descriptor f : CGAL::faces_around_target(halfedge(v, pmesh), pmesh)) + Vector_3 normal = cv_3(CGAL::NULL_VECTOR); + + halfedge_descriptor h = halfedge(v, pmesh); + if(h == boost::graph_traits::null_halfedge()) + return normal; + + halfedge_descriptor end = h; + do { - if(f == boost::graph_traits::null_face()) - continue; + if(!is_border(h, pmesh)) + { + if(vn_type == NO_WEIGHT) + { + const Vector_ref n = get(face_normals, face(h, pmesh)); + normal = traits.construct_sum_of_vectors_3_object()(normal, n); + } + else if(vn_type == SIN_WEIGHT) + { + const Vector_3 v1 = cv_3(get(vpmap, v), get(vpmap, source(h, pmesh))); + const Vector_3 v2 = cv_3(get(vpmap, v), get(vpmap, target(next(h, pmesh), pmesh))); - const Vector_ref n = get(face_normals, f); - normal = traits.construct_sum_of_vectors_3_object()(normal, n); + //v(i) and v(i+1) must me seen in ccw order, from v, so we reverse v1 and v2 + Vector_3 n = traits.construct_cross_product_vector_3_object()(v2, v1); + n = traits.construct_scaled_vector_3_object()(n, CGAL::approximate_sqrt(FT(1)/(csl_3(v1) * csl_3(v2)))); + + normal = traits.construct_sum_of_vectors_3_object()(normal, n); + } + else + { + std::cerr << "Error: unknown vertex normal type" << std::endl; + CGAL_assertion(false); + return CGAL::NULL_VECTOR; + } + } + + h = opposite(next(h, pmesh), pmesh); } + while(h != end); return normal; } @@ -609,7 +649,7 @@ compute_vertex_normal(typename boost::graph_traits::vertex_descript #ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP std::cout << "Failed to find most visible normal, use weighted sum of normals" << std::endl; #endif - normal = internal::compute_vertex_normal_as_weighted_sum_of_face_normals(v, face_normals, pmesh, traits); + normal = internal::compute_vertex_normal_as_sum_of_weighted_normals(v, internal::SIN_WEIGHT, face_normals, vpmap, pmesh, traits); } if(!traits.equal_3_object()(normal, CGAL::NULL_VECTOR)) From 6da2cebc8d01bc7d0f64cf44e4d619b4b03227f1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Mon, 26 Aug 2019 13:04:52 +0200 Subject: [PATCH 16/39] Sqrt -> Approximate_sqrt --- .../include/CGAL/Polygon_mesh_processing/compute_normal.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h index f73bb23492f..3e528a44408 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h @@ -259,7 +259,7 @@ bool does_enclose_other_normals(const int i, const int j, const int k, typename K::Compute_scalar_product_3 sp = traits.compute_scalar_product_3_object(); - const FT nbn = CGAL::sqrt(traits.compute_squared_length_3_object()(nb)); + const FT nbn = CGAL::approximate_sqrt(traits.compute_squared_length_3_object()(nb)); // check that this min circle defined by the diameter contains the other points const std::size_t nif = incident_faces.size(); From ae60ff2d15008cf5fab20eb28f7d28fa46d1c314 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Mon, 26 Aug 2019 13:05:27 +0200 Subject: [PATCH 17/39] Robustification of most visible vertex normal computation --- .../Polygon_mesh_processing/compute_normal.h | 54 +++++++++---------- 1 file changed, 27 insertions(+), 27 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h index 3e528a44408..3e454068c90 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h @@ -297,17 +297,21 @@ typename GT::Vector_3 compute_normals_bisector(const typename GT::Vector_3& ni, return traits.construct_sum_of_vectors_3_object()(ni, nj); // not normalized } -// @FIXME Handle the case of ni, nj and nk giving 3 points on a great circle template typename GT::Vector_3 compute_normals_bisector(const typename GT::Vector_3& ni, const typename GT::Vector_3& nj, const typename GT::Vector_3& nk, const GT& traits) { + typedef typename GT::FT FT; typedef typename GT::Point_3 Point_3; typedef typename GT::Vector_3 Vector_3; - Vector_3 nb = traits.construct_vector_3_object()(CGAL::NULL_VECTOR); + typename GT::Construct_scaled_vector_3 cslv_3 = traits.construct_scaled_vector_3_object(); + typename GT::Construct_sum_of_vectors_3 csv_3 = traits.construct_sum_of_vectors_3_object(); + typename GT::Construct_vector_3 cv_3 = traits.construct_vector_3_object(); + + Vector_3 nb = cv_3(CGAL::NULL_VECTOR); if(almost_equal(ni, nj, traits)) { @@ -329,27 +333,28 @@ typename GT::Vector_3 compute_normals_bisector(const typename GT::Vector_3& ni, CGAL_assertion(ni != nj); CGAL_assertion(ni != nk); CGAL_assertion(nj != nk); - CGAL_assertion(!traits.collinear_3_object()(CGAL::ORIGIN + ni, - CGAL::ORIGIN + nj, - CGAL::ORIGIN + nk)); + CGAL_assertion(!traits.collinear_3_object()(CGAL::ORIGIN + ni, CGAL::ORIGIN + nj, CGAL::ORIGIN + nk)); #ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP std::cout << "Triplet: ni[" << ni << "] nj[" << nj << "] nk[" << nk << "]" << std::endl; #endif - const Point_3 c = traits.construct_circumcenter_3_object()(CGAL::ORIGIN + ni, - CGAL::ORIGIN + nj, - CGAL::ORIGIN + nk); - CGAL_assertion(c != CGAL::ORIGIN); + const Point_3 c = traits.construct_circumcenter_3_object()(CGAL::ORIGIN + ni, CGAL::ORIGIN + nj, CGAL::ORIGIN + nk); + if(c == CGAL::ORIGIN) + { + // will happen if the three vectors live in the same plan, return some weighted sum + const FT third = FT(1)/FT(3); + return csv_3(csv_3(cslv_3(ni, third), cslv_3(nj, third)), cslv_3(nk, third)); + } - nb = traits.construct_vector_3_object()(CGAL::ORIGIN, c); // note that this isn't normalized + nb = cv_3(CGAL::ORIGIN, c); // note that this isn't normalized } return nb; } template -std::pair +typename GT::Vector_3 compute_most_visible_normal_2_points(std::vector::face_descriptor>& incident_faces, const FaceNormalVector& face_normals, const GT& traits) @@ -358,14 +363,15 @@ compute_most_visible_normal_2_points(std::vector::reference Vector_ref; - typename GT::Compute_scalar_product_3 sp = traits.compute_scalar_product_3_object(); + typename GT::Compute_scalar_product_3 sp_3 = traits.compute_scalar_product_3_object(); + typename GT::Construct_vector_3 cv_3 = traits.construct_vector_3_object(); #ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP std::cout << "Trying to find enclosing normal with 2 normals" << std::endl; #endif FT min_sp = -1; - Vector_3 n = traits.construct_vector_3_object()(CGAL::NULL_VECTOR); + Vector_3 n = cv_3(CGAL::NULL_VECTOR); const std::size_t nif = incident_faces.size(); for(int i=0; i= 0); if(sp_bi <= min_sp) @@ -396,8 +401,7 @@ compute_most_visible_normal_2_points(std::vector @@ -410,9 +414,6 @@ compute_most_visible_normal_3_points(const std::vector::reference Vector_ref; - typename GT::Compute_scalar_product_3 sp = traits.compute_scalar_product_3_object(); - typename GT::Construct_sum_of_vectors_3 csv = traits.construct_sum_of_vectors_3_object(); - #ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP std::cout << "Trying to find enclosing normal with 3 normals" << std::endl; #endif @@ -437,7 +438,7 @@ compute_most_visible_normal_3_points(const std::vector res = compute_most_visible_normal_2_points(incident_faces, face_normals, traits); + Vector_3 res = compute_most_visible_normal_2_points(incident_faces, face_normals, traits); - if(!res.second || // won't be able to find a most visible normal - (res.second && res.first != CGAL::NULL_VECTOR)) // found a valid normal through 2 point min circle - return res.first; + if(res != CGAL::NULL_VECTOR) // found a valid normal through 2 point min circle + return res; CGAL_assertion(incident_faces.size() > 2); From 5cf6ec52619b02c6c169d4b65a5d1ae8517c52b2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Mon, 26 Aug 2019 13:05:50 +0200 Subject: [PATCH 18/39] Fix minor bug in 'is_default_parameter' usage --- .../Polygon_mesh_processing/compute_normal.h | 57 +++++++++++-------- 1 file changed, 34 insertions(+), 23 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h index 3e454068c90..72a31119418 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h @@ -175,7 +175,7 @@ typename GetGeomTraits::type::Vector_3 compute_face_normal(typename boost::graph_traits::face_descriptor f, const PolygonMesh& pmesh) { - return compute_face_normal(f, pmesh, CGAL::Polygon_mesh_processing::parameters::all_default()); + return compute_face_normal(f, pmesh, CGAL::parameters::all_default()); } /** @@ -607,19 +607,20 @@ compute_vertex_normal(typename boost::graph_traits::vertex_descript typedef typename GT::Vector_3 Vector_3; GT traits = choose_parameter(get_parameter(np, internal_np::geom_traits), GT()); - typedef typename GetVertexPointMap::const_type VPMap; + typedef typename GetVertexPointMap::const_type VPMap; VPMap vpmap = choose_parameter(get_parameter(np, internal_np::vertex_point), get_const_property_map(vertex_point, pmesh)); typedef std::map Face_vector_map; typedef boost::associative_property_map Default_map; - typedef typename internal_np::Lookup_named_param_def ::type Face_normal_map; + typedef typename internal_np::Lookup_named_param_def::type Face_normal_map; + Face_vector_map default_fvmap; Face_normal_map face_normals = choose_parameter(get_parameter(np, internal_np::face_normal), - Default_map()); - const bool fnmap_valid = is_default_parameter(get_parameter(np, internal_np::face_normal)); + Default_map(default_fvmap)); + const bool must_compute_face_normals = is_default_parameter(get_parameter(np, internal_np::face_normal)); #ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP std::cout << std::endl << std::endl; @@ -632,7 +633,7 @@ compute_vertex_normal(typename boost::graph_traits::vertex_descript if(he == boost::graph_traits::null_halfedge()) return CGAL::NULL_VECTOR; - if(!fnmap_valid) + if(must_compute_face_normals) // if vertex normal type is sin-based weights, we don't need to compute normals { for(face_descriptor f : CGAL::faces_around_target(halfedge(v, pmesh), pmesh)) { @@ -663,16 +664,17 @@ typename GetGeomTraits::type::Vector_3 compute_vertex_normal(typename boost::graph_traits::vertex_descriptor v, const PolygonMesh& pmesh) { - return compute_vertex_normal(v, pmesh, CGAL::Polygon_mesh_processing::parameters::all_default()); + return compute_vertex_normal(v, pmesh, CGAL::parameters::all_default()); } /** * \ingroup PMP_normal_grp * computes the outward unit vector normal for all vertices of the polygon mesh. +* * @tparam PolygonMesh a model of `FaceListGraph` * @tparam VertexNormalMap a model of `WritablePropertyMap` with - `boost::graph_traits::%vertex_descriptor` as key type and - the return type of `compute_vertex_normal()` as value type. +* `boost::graph_traits::%vertex_descriptor` as key type and +* the return type of `compute_vertex_normal()` as value type. * * @param pmesh the polygon mesh * @param vnm the property map in which the normals are written @@ -698,22 +700,30 @@ void compute_vertex_normals(const PolygonMesh& pmesh, using parameters::is_default_parameter; using parameters::get_parameter; - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - typedef typename GetGeomTraits::type GT; - typedef typename GT::Vector_3 Vector_3; + typedef typename GetGeomTraits::type GT; + typedef typename GT::Vector_3 Vector_3; - typedef CGAL::dynamic_face_property_t Face_normal_tag; - typedef typename boost::property_map::type Face_normal_dmap; + typedef CGAL::dynamic_face_property_t Face_normal_tag; + typedef typename boost::property_map::const_type Face_normal_dmap; + +#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP + GT traits = choose_parameter(get_parameter(np, internal_np::geom_traits), GT()); + + typedef typename GetVertexPointMap::const_type VPMap; + VPMap vpmap = choose_parameter(get_parameter(np, internal_np::vertex_point), + get_const_property_map(vertex_point, pmesh)); +#endif typedef typename internal_np::Lookup_named_param_def::type Face_normal_map; + Face_normal_dmap>::type Face_normal_map; Face_normal_map face_normals = choose_parameter(get_parameter(np, internal_np::face_normal), get(Face_normal_tag(), pmesh)); - const bool fnmap_valid = is_default_parameter(get_parameter(np, internal_np::face_normal)); + const bool must_compute_face_normals = is_default_parameter(get_parameter(np, internal_np::face_normal)); - if(!fnmap_valid) + if(must_compute_face_normals) compute_face_normals(pmesh, face_normals, np); #ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG @@ -745,13 +755,14 @@ void compute_vertex_normals(const PolygonMesh& pmesh, VertexNormalMap vertex_nor /** * \ingroup PMP_normal_grp * computes the outward unit vector normal for all vertices and faces of the polygon mesh. +* * @tparam PolygonMesh a model of `FaceListGraph` * @tparam VertexNormalMap a model of `WritablePropertyMap` with - `boost::graph_traits::%vertex_descriptor` as key type and - `Kernel::Vector_3` as value type. +* `boost::graph_traits::%vertex_descriptor` as key type and +* `Kernel::Vector_3` as value type. * @tparam FaceNormalMap a model of `ReadWritePropertyMap` with - `boost::graph_traits::%face_descriptor` as key type and - `Kernel::Vector_3` as value type. +* `boost::graph_traits::%face_descriptor` as key type and +* `Kernel::Vector_3` as value type. * * @param pmesh the polygon mesh * @param vertex_normals the property map in which the vertex normals are written From d9c73a1b35a9d0d65f0d7d1753cf7d32bd2a31f6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Mon, 26 Aug 2019 13:06:15 +0200 Subject: [PATCH 19/39] Properly test the family of normal computation functions --- .../data/folded_star.off | 15 ++ .../pmp_compute_normals_test.cpp | 166 ++++++++++++++---- 2 files changed, 142 insertions(+), 39 deletions(-) create mode 100644 Polygon_mesh_processing/test/Polygon_mesh_processing/data/folded_star.off diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/data/folded_star.off b/Polygon_mesh_processing/test/Polygon_mesh_processing/data/folded_star.off new file mode 100644 index 00000000000..68c8a887da2 --- /dev/null +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/data/folded_star.off @@ -0,0 +1,15 @@ +OFF +7 6 0 +0.50085999554530469 0.98935974982034269 0.054777016186685679 +0.51086414037695982 0.0074532224578913869 0.0098759303387234813 +0.23667027533904861 -0.19236191791372245 0.078118737677961653 +0.49890739640011506 0.47147692553082188 0.51901846156493603 +0.5 0.5 0 +0.8400145559724278 0.28990007922508693 -0.53104077094013602 +-0.011501240589425327 0.99837637926383715 0.011228717925938695 +3 4 5 0 +3 1 4 3 +3 3 4 6 +3 4 0 6 +3 1 2 4 +3 2 5 4 diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp index 364df9d3cfa..4d222f79c3c 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp @@ -1,62 +1,150 @@ #include #include -#include -#include -#include +#include +#include + +#include #include #include -typedef CGAL::Exact_predicates_inexact_constructions_kernel Epic; -typedef CGAL::Exact_predicates_exact_constructions_kernel Epec; +typedef CGAL::Exact_predicates_inexact_constructions_kernel EPICK; +typedef CGAL::Exact_predicates_exact_constructions_kernel EPECK; -template -void test_normals(const char* file_name) +typedef CGAL::Surface_mesh EPICK_SM; +typedef CGAL::Surface_mesh EPECK_SM; + +namespace PMP = CGAL::Polygon_mesh_processing; + +template +void test(const Mesh& mesh, + const VertexNormalPmap& vnormals, + const FaceNormalPmap& fnormals) { - typedef typename K::Point_3 Point; - typedef typename K::Vector_3 Vector; - typedef CGAL::Surface_mesh Surface_mesh; - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - typedef typename boost::graph_traits::face_descriptor face_descriptor; - - Surface_mesh mesh; + typedef typename K::Vector_3 Vector; + + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::face_descriptor face_descriptor; + + typedef typename PMP::GetVertexPointMap::const_type VPMap; + VPMap vpmap = get_const_property_map(CGAL::vertex_point, mesh); + + const vertex_descriptor first_vertex = *(vertices(mesh).begin()); + const face_descriptor first_face = *(faces(mesh).begin()); + + PMP::compute_face_normals(mesh, fnormals); + PMP::compute_face_normals(mesh, fnormals, PMP::parameters::vertex_point_map(vpmap)); + PMP::compute_face_normals(mesh, fnormals, CGAL::parameters::vertex_point_map(vpmap) + .geom_traits(K())); + + Vector f0n = PMP::compute_face_normal(first_face, mesh); + assert(f0n == get(fnormals, first_face)); + PMP::compute_face_normal(first_face, mesh, PMP::parameters::vertex_point_map(vpmap)); + + PMP::compute_vertex_normals(mesh, vnormals); + PMP::compute_vertex_normals(mesh, vnormals, PMP::parameters::vertex_point_map(vpmap)); + PMP::compute_vertex_normals(mesh, vnormals, CGAL::parameters::vertex_point_map(vpmap) + .geom_traits(K())); + + Vector v0n = PMP::compute_vertex_normal(first_vertex, mesh); + assert(v0n == get(vnormals, first_vertex)); + v0n = PMP::compute_vertex_normal(first_vertex, mesh, PMP::parameters::vertex_point_map(vpmap) + .face_normal_map(fnormals)); + assert(v0n == get(vnormals, first_vertex)); + + PMP::compute_normals(mesh, vnormals, fnormals); + PMP::compute_normals(mesh, vnormals, fnormals, PMP::parameters::vertex_point_map(vpmap)); + PMP::compute_normals(mesh, vnormals, fnormals, CGAL::parameters::vertex_point_map(vpmap) + .geom_traits(K())); + + for(vertex_descriptor v : vertices(mesh)) { + assert(get(vnormals, v) != CGAL::NULL_VECTOR); + } + + for(face_descriptor f : faces(mesh)) { + assert(get(fnormals, f) != CGAL::NULL_VECTOR); + } +} + +template +void test_SM(const char* file_name) +{ + typedef CGAL::Surface_mesh SM; + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::face_descriptor face_descriptor; + + typedef typename K::Vector_3 Vector; + + SM mesh; std::ifstream input(file_name); - if (!(input >> mesh)){ - std::cerr << "Error: cannot read Surface_mesh : " << file_name << "\n"; + if(!(input >> mesh)) + { + std::cerr << "Error: cannot read " << file_name << " as a Surface_mesh\n"; assert(false); } - - typename Surface_mesh::template Property_map fnormals; - bool created; - boost::tie(fnormals, created) = mesh.template add_property_map("f:normals",Vector(0,0,0)); - CGAL::Polygon_mesh_processing::compute_face_normals(mesh, fnormals); - CGAL::Polygon_mesh_processing::compute_face_normals(mesh, fnormals, - CGAL::Polygon_mesh_processing::parameters::vertex_point_map(mesh.points())); - CGAL::Polygon_mesh_processing::compute_face_normals(mesh, fnormals, - CGAL::Polygon_mesh_processing::parameters::vertex_point_map(mesh.points()).geom_traits(K())); - typename Surface_mesh:: template Property_map vnormals; + typename SM::template Property_map vnormals; + vnormals = mesh.template add_property_map("v:normals", CGAL::NULL_VECTOR).first; + typename SM::template Property_map fnormals; + fnormals = mesh.template add_property_map("f:normals", CGAL::NULL_VECTOR).first; - boost::tie(vnormals, created) = mesh.template add_property_map("v:normals",Vector(0,0,0)); - CGAL::Polygon_mesh_processing::compute_vertex_normals(mesh, vnormals); - CGAL::Polygon_mesh_processing::compute_vertex_normals(mesh, vnormals, - CGAL::Polygon_mesh_processing::parameters::vertex_point_map(mesh.points())); - CGAL::Polygon_mesh_processing::compute_vertex_normals(mesh, vnormals, - CGAL::Polygon_mesh_processing::parameters::vertex_point_map(mesh.points()).geom_traits(K())); + test(mesh, vnormals, fnormals); +} - CGAL::Polygon_mesh_processing::compute_normals(mesh, vnormals, fnormals); - CGAL::Polygon_mesh_processing::compute_normals(mesh, vnormals, fnormals, - CGAL::Polygon_mesh_processing::parameters::vertex_point_map(mesh.points())); - CGAL::Polygon_mesh_processing::compute_normals(mesh, vnormals, fnormals, - CGAL::Polygon_mesh_processing::parameters::vertex_point_map(mesh.points()).geom_traits(K())); +template +void test_Polyhedron(const char* file_name) +{ + typedef CGAL::Polyhedron_3 Polyhedron; + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::face_descriptor face_descriptor; + typedef typename K::Vector_3 Vector; + + Polyhedron mesh; + std::ifstream input(file_name); + if(!(input >> mesh)) + { + std::cerr << "Error: cannot read " << file_name << " as a Polyhedron\n"; + assert(false); + } + + typedef std::map VertexNormalMap; + typedef boost::associative_property_map VertexNormalPMap; + typedef std::map FaceNormalMap; + typedef boost::associative_property_map FaceNormalPMap; + + VertexNormalMap vn_map; + FaceNormalMap fn_map; + + for(vertex_descriptor v : vertices(mesh)) + vn_map[v] = CGAL::NULL_VECTOR; + for(face_descriptor f : faces(mesh)) + fn_map[f] = CGAL::NULL_VECTOR; + + VertexNormalPMap vnormals(vn_map); + FaceNormalPMap fnormals(fn_map); + + test(mesh, vnormals, fnormals); +} + +void test(const char* filename) +{ + std::cout << "test " << filename << "..." << std::endl; + + test_SM(filename); + test_SM(filename); + test_Polyhedron(filename); + test_Polyhedron(filename); } int main() { - test_normals("data/elephant.off"); - test_normals("data/elephant.off"); + test("../data/elephant.off"); + test("../data/folded_star.off"); + test("../data/joint_refined.off"); + test("../data/mannequin-devil.off"); + test("../data/U.off"); std::cerr << "All done." << std::endl; return EXIT_SUCCESS; From 9ea2f668a6d1581188a219d9af3de5fa118cf346 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Mon, 26 Aug 2019 13:19:40 +0200 Subject: [PATCH 20/39] Fix data file links --- .../pmp_compute_normals_test.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp index 4d222f79c3c..c1942339d09 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp @@ -140,11 +140,11 @@ void test(const char* filename) int main() { - test("../data/elephant.off"); - test("../data/folded_star.off"); - test("../data/joint_refined.off"); - test("../data/mannequin-devil.off"); - test("../data/U.off"); + test("data/elephant.off"); + test("data/folded_star.off"); + test("data/joint_refined.off"); + test("data/mannequin-devil.off"); + test("data/U.off"); std::cerr << "All done." << std::endl; return EXIT_SUCCESS; From 61ef88f2de48b67ed56904c86347f96f5c1a0c20 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Mon, 26 Aug 2019 13:24:43 +0200 Subject: [PATCH 21/39] Clean example --- .../compute_normals_example.cpp | 39 +++++++++---------- 1 file changed, 18 insertions(+), 21 deletions(-) diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/compute_normals_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/compute_normals_example.cpp index ac4e2c3b4a3..cd7ba66f2e4 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/compute_normals_example.cpp +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/compute_normals_example.cpp @@ -6,13 +6,16 @@ #include #include -typedef CGAL::Exact_predicates_inexact_constructions_kernel K; -typedef K::Point_3 Point; -typedef K::Vector_3 Vector; +namespace PMP = CGAL::Polygon_mesh_processing; -typedef CGAL::Surface_mesh Surface_mesh; -typedef boost::graph_traits::vertex_descriptor vertex_descriptor; -typedef boost::graph_traits::face_descriptor face_descriptor; +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; + +typedef K::Point_3 Point; +typedef K::Vector_3 Vector; + +typedef CGAL::Surface_mesh Surface_mesh; +typedef boost::graph_traits::vertex_descriptor vertex_descriptor; +typedef boost::graph_traits::face_descriptor face_descriptor; int main(int argc, char* argv[]) { @@ -25,24 +28,18 @@ int main(int argc, char* argv[]) return 1; } - auto fnormals = mesh.add_property_map - ("f:normals", CGAL::NULL_VECTOR).first; - auto vnormals = mesh.add_property_map - ("v:normals", CGAL::NULL_VECTOR).first; + auto vnormals = mesh.add_property_map("v:normals", CGAL::NULL_VECTOR).first; + auto fnormals = mesh.add_property_map("f:normals", CGAL::NULL_VECTOR).first; - CGAL::Polygon_mesh_processing::compute_normals(mesh, - vnormals, - fnormals, - CGAL::Polygon_mesh_processing::parameters::vertex_point_map(mesh.points()). - geom_traits(K())); + PMP::compute_normals(mesh, vnormals, fnormals); + + std::cout << "Vertex normals :" << std::endl; + for(vertex_descriptor vd: vertices(mesh)) + std::cout << vnormals[vd] << std::endl; std::cout << "Face normals :" << std::endl; - for(face_descriptor fd: faces(mesh)){ + for(face_descriptor fd: faces(mesh)) std::cout << fnormals[fd] << std::endl; - } - std::cout << "Vertex normals :" << std::endl; - for(vertex_descriptor vd: vertices(mesh)){ - std::cout << vnormals[vd] << std::endl; - } + return 0; } From c2cb44175b73b32b1573c01e405dba98a8461543 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Mon, 26 Aug 2019 13:24:51 +0200 Subject: [PATCH 22/39] Add a reference to the paper --- .../include/CGAL/Polygon_mesh_processing/compute_normal.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h index 72a31119418..ac0351b6811 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h @@ -464,7 +464,7 @@ compute_most_visible_normal_3_points(const std::vector typename GT::Vector_3 compute_vertex_normal_most_visible_min_circle(typename boost::graph_traits::vertex_descriptor v, From 2c06609445d23a19275f269956ab6bd4a295947b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Fri, 30 Aug 2019 16:22:50 +0200 Subject: [PATCH 23/39] Fix some warnings --- .../Polygon_mesh_processing/compute_normal.h | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h index ac0351b6811..dacf2b532b3 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h @@ -246,7 +246,7 @@ bool almost_equal(const typename GT::Vector_3& v1, const typename GT::Vector_3& } template -bool does_enclose_other_normals(const int i, const int j, const int k, +bool does_enclose_other_normals(const std::size_t i, const std::size_t j, const std::size_t k, const typename K::Vector_3& nb, const typename K::FT sp_bi, const std::vector::face_descriptor>& incident_faces, @@ -254,7 +254,6 @@ bool does_enclose_other_normals(const int i, const int j, const int k, const K& traits) { typedef typename K::FT FT; - typedef typename K::Vector_3 Vector_3; typedef typename boost::property_traits::reference Vector_ref; typename K::Compute_scalar_product_3 sp = traits.compute_scalar_product_3_object(); @@ -263,7 +262,7 @@ bool does_enclose_other_normals(const int i, const int j, const int k, // check that this min circle defined by the diameter contains the other points const std::size_t nif = incident_faces.size(); - for(int l=0; l::face_descriptor face_descriptor; - typedef typename GT::FT FT; typedef typename GT::Vector_3 Vector_3; std::vector incident_faces; From dff631f2bffabf0a4ee12c5f32f59d87d87b0370 Mon Sep 17 00:00:00 2001 From: Mael Date: Fri, 15 Nov 2019 10:14:21 +0100 Subject: [PATCH 24/39] Fix namespace --- .../test/Polygon_mesh_processing/pmp_compute_normals_test.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp index c1942339d09..2485a85ca93 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp @@ -27,7 +27,7 @@ void test(const Mesh& mesh, typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; typedef typename boost::graph_traits::face_descriptor face_descriptor; - typedef typename PMP::GetVertexPointMap::const_type VPMap; + typedef typename CGAL::GetVertexPointMap::const_type VPMap; VPMap vpmap = get_const_property_map(CGAL::vertex_point, mesh); const vertex_descriptor first_vertex = *(vertices(mesh).begin()); From a47f0ee91d7742ae411f2e94f76041009e560236 Mon Sep 17 00:00:00 2001 From: Mael Date: Fri, 15 Nov 2019 12:03:18 +0100 Subject: [PATCH 25/39] Fix inconsistency in doc/code parameter names --- .../include/CGAL/Polygon_mesh_processing/compute_normal.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h index 237b6e98095..f7558718e74 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h @@ -666,7 +666,7 @@ compute_vertex_normal(typename boost::graph_traits::vertex_descript * the return type of `compute_vertex_normal()` as value type. * * @param pmesh the polygon mesh -* @param vnm the property map in which the normals are written +* @param vertex_normals the property map in which the normals are written * @param np optional sequence of \ref pmp_namedparameters "Named Parameters" among the ones listed below * * \cgalNamedParamsBegin From 0ddece5c81b2f4f2079319f3f18275ca5e04c8cf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Mon, 25 Nov 2019 09:43:41 +0100 Subject: [PATCH 26/39] Try to fix some weird assertion in a 32 bit test --- .../test/Polygon_mesh_processing/pmp_compute_normals_test.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp index 2485a85ca93..101e2c6fcca 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp @@ -140,6 +140,8 @@ void test(const char* filename) int main() { + CGAL::Set_ieee_double_precision pfr; + test("data/elephant.off"); test("data/folded_star.off"); test("data/joint_refined.off"); From e50b54a8b8c1410904492d19795de757f98c71d7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Thu, 28 Nov 2019 18:28:19 +0100 Subject: [PATCH 27/39] Minor cleaning --- .../include/CGAL/Polygon_mesh_processing/compute_normal.h | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h index f7558718e74..1fcc36ce51b 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h @@ -614,7 +614,8 @@ compute_vertex_normal(typename boost::graph_traits::vertex_descript #ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP std::cout << std::endl << std::endl; std::cout << "----------------------------------------------------------------------" << std::endl; - std::cout << "compute vertex at " << get(vpmap, v) << std::endl; + std::cout << "compute vertex at " << get(vpmap, v) + << ", must compute face normals? " << must_compute_face_normals << std::endl; #endif // handle isolated vertices @@ -622,7 +623,7 @@ compute_vertex_normal(typename boost::graph_traits::vertex_descript if(he == boost::graph_traits::null_halfedge()) return CGAL::NULL_VECTOR; - if(must_compute_face_normals) // if vertex normal type is sin-based weights, we don't need to compute normals + if(must_compute_face_normals) { for(face_descriptor f : CGAL::faces_around_target(halfedge(v, pmesh), pmesh)) { @@ -639,7 +640,8 @@ compute_vertex_normal(typename boost::graph_traits::vertex_descript #ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP std::cout << "Failed to find most visible normal, use weighted sum of normals" << std::endl; #endif - normal = internal::compute_vertex_normal_as_sum_of_weighted_normals(v, internal::SIN_WEIGHT, face_normals, vpmap, pmesh, traits); + normal = internal::compute_vertex_normal_as_sum_of_weighted_normals( + v, internal::SIN_WEIGHT, face_normals, vpmap, pmesh, traits); } if(!traits.equal_3_object()(normal, CGAL::NULL_VECTOR)) From 908da8cdf78538e4f022567a2aee0c0b5494d2f2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Thu, 28 Nov 2019 18:39:12 +0100 Subject: [PATCH 28/39] Fix some debug macros --- .../include/CGAL/Polygon_mesh_processing/compute_normal.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h index 1fcc36ce51b..46e12778ce1 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h @@ -699,7 +699,7 @@ void compute_vertex_normals(const PolygonMesh& pmesh, typedef CGAL::dynamic_face_property_t Face_normal_tag; typedef typename boost::property_map::const_type Face_normal_dmap; -#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP +#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG GT traits = choose_parameter(get_parameter(np, internal_np::geom_traits), GT()); typedef typename GetVertexPointMap::const_type VPMap; @@ -717,7 +717,7 @@ void compute_vertex_normals(const PolygonMesh& pmesh, if(must_compute_face_normals) compute_face_normals(pmesh, face_normals, np); -#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG +#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP std::ofstream out("computed_normals.cgal.polylines.txt"); const Bbox_3 bb = bbox(pmesh, np); const typename GT::FT bbox_diagonal = CGAL::sqrt(CGAL::square(bb.xmax() - bb.xmin()) + @@ -730,7 +730,7 @@ void compute_vertex_normals(const PolygonMesh& pmesh, const Vector_3 n = compute_vertex_normal(v, pmesh, np.face_normal_map(face_normals)); put(vertex_normals, v, n); -#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG +#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP out << "2 " << get(vpmap, v) << " " << get(vpmap, v) + traits.construct_scaled_vector_3_object()(n, 0.1 * bbox_diagonal) << "\n"; #endif From 783d5df98cc6c9457e05f034051a501b14a34a85 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Thu, 28 Nov 2019 18:39:28 +0100 Subject: [PATCH 29/39] TMP DEBUG CODE --- .../include/CGAL/Polygon_mesh_processing/compute_normal.h | 2 +- .../test/Polygon_mesh_processing/pmp_compute_normals_test.cpp | 4 ++++ 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h index 46e12778ce1..76687bd1776 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h @@ -611,7 +611,7 @@ compute_vertex_normal(typename boost::graph_traits::vertex_descript Default_map(default_fvmap)); const bool must_compute_face_normals = is_default_parameter(get_parameter(np, internal_np::face_normal)); -#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP +#if 1//def CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP std::cout << std::endl << std::endl; std::cout << "----------------------------------------------------------------------" << std::endl; std::cout << "compute vertex at " << get(vpmap, v) diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp index 101e2c6fcca..cbe56ae2b75 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp @@ -1,3 +1,5 @@ +#define CGAL_PMP_COMPUTE_NORMAL_DEBUG + #include #include @@ -51,6 +53,8 @@ void test(const Mesh& mesh, assert(v0n == get(vnormals, first_vertex)); v0n = PMP::compute_vertex_normal(first_vertex, mesh, PMP::parameters::vertex_point_map(vpmap) .face_normal_map(fnormals)); + std::cout.precision(17); + std::cout << v0n << " versus " << get(vnormals, first_vertex) << std::endl; assert(v0n == get(vnormals, first_vertex)); PMP::compute_normals(mesh, vnormals, fnormals); From b90128af9dab9f61f0fe1bf27c76d7507ea32617 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Fri, 6 Dec 2019 09:14:27 +0100 Subject: [PATCH 30/39] Fix unused warning --- .../include/CGAL/Polygon_mesh_processing/compute_normal.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h index 76687bd1776..681ed0b4c64 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h @@ -699,7 +699,7 @@ void compute_vertex_normals(const PolygonMesh& pmesh, typedef CGAL::dynamic_face_property_t Face_normal_tag; typedef typename boost::property_map::const_type Face_normal_dmap; -#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG +#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP GT traits = choose_parameter(get_parameter(np, internal_np::geom_traits), GT()); typedef typename GetVertexPointMap::const_type VPMap; From 257ad60fef14aefeb87c7787ee0d761caa01c644 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Tue, 10 Dec 2019 16:36:16 +0100 Subject: [PATCH 31/39] More debug code to understand Debian 32 --- .../CGAL/Polygon_mesh_processing/compute_normal.h | 13 +++++++++++++ .../pmp_compute_normals_test.cpp | 14 ++++++++++---- 2 files changed, 23 insertions(+), 4 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h index 681ed0b4c64..c027795a773 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h @@ -87,13 +87,19 @@ void sum_normals(const PM& pmesh, typedef typename boost::property_traits::reference Point_ref; + std::cout << "sum_normals..." << std::endl; + halfedge_descriptor he = halfedge(f, pmesh); vertex_descriptor v = source(he, pmesh); const Point_ref pv = get(vpmap, v); + std::cout << "pv: " << pv << std::endl; + while(v != target(next(he, pmesh), pmesh)) { const Point_ref pvn = get(vpmap, target(he, pmesh)); const Point_ref pvnn = get(vpmap, target(next(he, pmesh), pmesh)); + std::cout << "pvn: " << pvn << std::endl; + std::cout << "pvnn: " << pvnn << std::endl; const Vector n = internal::triangle_normal(pv, pvn, pvnn, traits); sum = traits.construct_sum_of_vectors_3_object()(sum, n); @@ -634,6 +640,13 @@ compute_vertex_normal(typename boost::graph_traits::vertex_descript } } + std::cout << "Incident face normals:" << std::endl; + for(halfedge_descriptor h : CGAL::halfedges_around_target(v, pmesh)) + { + if(!is_border(h, pmesh)) + std::cout << "normal: " << get(face_normals, face(h, pmesh)) << std::endl; + } + Vector_3 normal = internal::compute_vertex_normal_most_visible_min_circle(v, face_normals, pmesh, traits); if(traits.equal_3_object()(normal, CGAL::NULL_VECTOR)) // can't always find a most visible normal { diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp index cbe56ae2b75..9b61436b6e2 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp @@ -32,6 +32,8 @@ void test(const Mesh& mesh, typedef typename CGAL::GetVertexPointMap::const_type VPMap; VPMap vpmap = get_const_property_map(CGAL::vertex_point, mesh); + assert(!is_empty(mesh)); + const vertex_descriptor first_vertex = *(vertices(mesh).begin()); const face_descriptor first_face = *(faces(mesh).begin()); @@ -80,6 +82,8 @@ void test_SM(const char* file_name) typedef typename K::Vector_3 Vector; + std::cout << "Test with Surface_mesh, and kernel: " << typeid(K).name() << std::endl; + SM mesh; std::ifstream input(file_name); if(!(input >> mesh)) @@ -105,6 +109,8 @@ void test_Polyhedron(const char* file_name) typedef typename K::Vector_3 Vector; + std::cout << "Test with Polyhedron, and kernel: " << typeid(K).name() << std::endl; + Polyhedron mesh; std::ifstream input(file_name); if(!(input >> mesh)) @@ -146,10 +152,10 @@ int main() { CGAL::Set_ieee_double_precision pfr; - test("data/elephant.off"); - test("data/folded_star.off"); - test("data/joint_refined.off"); - test("data/mannequin-devil.off"); +// test("data/elephant.off"); +// test("data/folded_star.off"); +// test("data/joint_refined.off"); +// test("data/mannequin-devil.off"); test("data/U.off"); std::cerr << "All done." << std::endl; From 1a4b8e8dc0b2aeb006061abc831bd4e787d40edb Mon Sep 17 00:00:00 2001 From: Mael Date: Wed, 18 Dec 2019 09:09:38 +0100 Subject: [PATCH 32/39] more debug code for Debian 32 --- .../include/CGAL/Polygon_mesh_processing/compute_normal.h | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h index c027795a773..9e62c0ff12b 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h @@ -159,12 +159,16 @@ compute_face_normal(typename boost::graph_traits::face_descriptor f typedef typename GT::Point_3 Point; typedef typename GT::Vector_3 Vector_3; + std::cout << "Call to compute_face_normal()" << std::endl; + Vector_3 normal = traits.construct_vector_3_object()(CGAL::NULL_VECTOR); internal::sum_normals(pmesh, f, vpmap, normal, traits); if(!traits.equal_3_object()(normal, CGAL::NULL_VECTOR)) internal::normalize(normal, traits); + std::cout << "result: " << normal << std::endl; + return normal; } From 79756e56cd9b660a18ac504173fa0f858f8e8110 Mon Sep 17 00:00:00 2001 From: Mael Date: Wed, 18 Dec 2019 09:10:43 +0100 Subject: [PATCH 33/39] Disable tests that work (tmp - will be removed) --- .../Polygon_mesh_processing/pmp_compute_normals_test.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp index 9b61436b6e2..b7c6c557d3b 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp @@ -142,10 +142,10 @@ void test(const char* filename) { std::cout << "test " << filename << "..." << std::endl; - test_SM(filename); +// test_SM(filename); test_SM(filename); - test_Polyhedron(filename); - test_Polyhedron(filename); +// test_Polyhedron(filename); +// test_Polyhedron(filename); } int main() From b005350f375aaffd8861ad500b000aaf07962043 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Mon, 30 Dec 2019 11:25:21 +0100 Subject: [PATCH 34/39] more debug code, exact 32bit doesn't seem deterministic --- .../include/CGAL/Polygon_mesh_processing/compute_normal.h | 5 ++++- .../Polygon_mesh_processing/pmp_compute_normals_test.cpp | 2 ++ 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h index 9e62c0ff12b..32c11c2cbd5 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h @@ -214,6 +214,7 @@ void compute_face_normals(const PolygonMesh& pmesh, { typename Kernel::Vector_3 vec = compute_face_normal(f, pmesh, np); put(face_normals, f, vec); + std::cout << "normal at face " << f << " is " << get(face_normals, f) << std::endl; } } @@ -644,11 +645,13 @@ compute_vertex_normal(typename boost::graph_traits::vertex_descript } } + std::cout << "type of FT: " << typeid(typename GT::FT).name() << std::endl; + std::cout << "Incident face normals:" << std::endl; for(halfedge_descriptor h : CGAL::halfedges_around_target(v, pmesh)) { if(!is_border(h, pmesh)) - std::cout << "normal: " << get(face_normals, face(h, pmesh)) << std::endl; + std::cout << "get normal at f " << face(h, pmesh) << " : " << get(face_normals, face(h, pmesh)).approx() << std::endl; } Vector_3 normal = internal::compute_vertex_normal_most_visible_min_circle(v, face_normals, pmesh, traits); diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp index b7c6c557d3b..80200c73ddc 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp @@ -150,6 +150,8 @@ void test(const char* filename) int main() { + std::cout.precision(17); + CGAL::Set_ieee_double_precision pfr; // test("data/elephant.off"); From a53a66013d7e77331a001839af061d23d3f5e2df Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Thu, 2 Jan 2020 09:32:43 +0100 Subject: [PATCH 35/39] Add some verbose macro --- .../include/CGAL/Polygon_mesh_processing/compute_normal.h | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h index 32c11c2cbd5..cb7c4bd8390 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h @@ -214,7 +214,9 @@ void compute_face_normals(const PolygonMesh& pmesh, { typename Kernel::Vector_3 vec = compute_face_normal(f, pmesh, np); put(face_normals, f, vec); +#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG std::cout << "normal at face " << f << " is " << get(face_normals, f) << std::endl; +#endif } } @@ -647,12 +649,14 @@ compute_vertex_normal(typename boost::graph_traits::vertex_descript std::cout << "type of FT: " << typeid(typename GT::FT).name() << std::endl; +#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG std::cout << "Incident face normals:" << std::endl; for(halfedge_descriptor h : CGAL::halfedges_around_target(v, pmesh)) { if(!is_border(h, pmesh)) - std::cout << "get normal at f " << face(h, pmesh) << " : " << get(face_normals, face(h, pmesh)).approx() << std::endl; + std::cout << "get normal at f " << face(h, pmesh) << " : " << get(face_normals, face(h, pmesh)) << std::endl; } +#endif Vector_3 normal = internal::compute_vertex_normal_most_visible_min_circle(v, face_normals, pmesh, traits); if(traits.equal_3_object()(normal, CGAL::NULL_VECTOR)) // can't always find a most visible normal From 151375255902c6b954237c9ab2ef4edfcf7a2f28 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Thu, 2 Jan 2020 09:33:41 +0100 Subject: [PATCH 36/39] Switch the test to exact square root computations --- .../pmp_compute_normals_test.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp index 80200c73ddc..219742bdd99 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp @@ -1,7 +1,7 @@ #define CGAL_PMP_COMPUTE_NORMAL_DEBUG #include -#include +#include #include #include @@ -11,11 +11,11 @@ #include #include -typedef CGAL::Exact_predicates_inexact_constructions_kernel EPICK; -typedef CGAL::Exact_predicates_exact_constructions_kernel EPECK; +typedef CGAL::Exact_predicates_inexact_constructions_kernel EPICK; +typedef CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt EPECK; -typedef CGAL::Surface_mesh EPICK_SM; -typedef CGAL::Surface_mesh EPECK_SM; +typedef CGAL::Surface_mesh EPICK_SM; +typedef CGAL::Surface_mesh EPECK_SM; namespace PMP = CGAL::Polygon_mesh_processing; From e388e63e0b5d54c0c855d505343670e986587f88 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Mon, 6 Jan 2020 10:53:13 +0100 Subject: [PATCH 37/39] Require Core for PMP/test/pmp_compute_normals_test --- .../test/Polygon_mesh_processing/CMakeLists.txt | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt index 02f55be25dd..2bf92f35e85 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt @@ -5,7 +5,7 @@ cmake_minimum_required(VERSION 3.1...3.15) project( Polygon_mesh_processing_Tests ) # CGAL and its components -find_package( CGAL QUIET COMPONENTS ) +find_package( CGAL QUIET COMPONENTS Core) if ( NOT CGAL_FOUND ) @@ -14,6 +14,8 @@ if ( NOT CGAL_FOUND ) endif() +include(${CGAL_USE_FILE}) + # Boost and its components find_package( Boost REQUIRED ) From eeeb5929cae09c00d1b187f07a3a760c693815d2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Mon, 6 Jan 2020 10:53:51 +0100 Subject: [PATCH 38/39] Clean debug code --- .../Polygon_mesh_processing/compute_normal.h | 13 +------------ .../pmp_compute_normals_test.cpp | 17 ++++++++--------- 2 files changed, 9 insertions(+), 21 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h index cb7c4bd8390..e3696de26ff 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h @@ -87,19 +87,14 @@ void sum_normals(const PM& pmesh, typedef typename boost::property_traits::reference Point_ref; - std::cout << "sum_normals..." << std::endl; - halfedge_descriptor he = halfedge(f, pmesh); vertex_descriptor v = source(he, pmesh); const Point_ref pv = get(vpmap, v); - std::cout << "pv: " << pv << std::endl; while(v != target(next(he, pmesh), pmesh)) { const Point_ref pvn = get(vpmap, target(he, pmesh)); const Point_ref pvnn = get(vpmap, target(next(he, pmesh), pmesh)); - std::cout << "pvn: " << pvn << std::endl; - std::cout << "pvnn: " << pvnn << std::endl; const Vector n = internal::triangle_normal(pv, pvn, pvnn, traits); sum = traits.construct_sum_of_vectors_3_object()(sum, n); @@ -159,16 +154,12 @@ compute_face_normal(typename boost::graph_traits::face_descriptor f typedef typename GT::Point_3 Point; typedef typename GT::Vector_3 Vector_3; - std::cout << "Call to compute_face_normal()" << std::endl; - Vector_3 normal = traits.construct_vector_3_object()(CGAL::NULL_VECTOR); internal::sum_normals(pmesh, f, vpmap, normal, traits); if(!traits.equal_3_object()(normal, CGAL::NULL_VECTOR)) internal::normalize(normal, traits); - std::cout << "result: " << normal << std::endl; - return normal; } @@ -624,7 +615,7 @@ compute_vertex_normal(typename boost::graph_traits::vertex_descript Default_map(default_fvmap)); const bool must_compute_face_normals = is_default_parameter(get_parameter(np, internal_np::face_normal)); -#if 1//def CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP +#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP std::cout << std::endl << std::endl; std::cout << "----------------------------------------------------------------------" << std::endl; std::cout << "compute vertex at " << get(vpmap, v) @@ -647,8 +638,6 @@ compute_vertex_normal(typename boost::graph_traits::vertex_descript } } - std::cout << "type of FT: " << typeid(typename GT::FT).name() << std::endl; - #ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG std::cout << "Incident face normals:" << std::endl; for(halfedge_descriptor h : CGAL::halfedges_around_target(v, pmesh)) diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp index 219742bdd99..cd49763c1f7 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp @@ -1,4 +1,4 @@ -#define CGAL_PMP_COMPUTE_NORMAL_DEBUG +// #define CGAL_PMP_COMPUTE_NORMAL_DEBUG #include #include @@ -56,7 +56,6 @@ void test(const Mesh& mesh, v0n = PMP::compute_vertex_normal(first_vertex, mesh, PMP::parameters::vertex_point_map(vpmap) .face_normal_map(fnormals)); std::cout.precision(17); - std::cout << v0n << " versus " << get(vnormals, first_vertex) << std::endl; assert(v0n == get(vnormals, first_vertex)); PMP::compute_normals(mesh, vnormals, fnormals); @@ -142,10 +141,10 @@ void test(const char* filename) { std::cout << "test " << filename << "..." << std::endl; -// test_SM(filename); + test_SM(filename); test_SM(filename); -// test_Polyhedron(filename); -// test_Polyhedron(filename); + test_Polyhedron(filename); + test_Polyhedron(filename); } int main() @@ -154,10 +153,10 @@ int main() CGAL::Set_ieee_double_precision pfr; -// test("data/elephant.off"); -// test("data/folded_star.off"); -// test("data/joint_refined.off"); -// test("data/mannequin-devil.off"); + test("data/elephant.off"); + test("data/folded_star.off"); + test("data/joint_refined.off"); + test("data/mannequin-devil.off"); test("data/U.off"); std::cerr << "All done." << std::endl; From 6dce5f3115aa1e1a53c138a27e7465ceb1703aca Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Mon, 6 Jan 2020 11:05:00 +0100 Subject: [PATCH 39/39] Disable tests using EPECK (for performance reasons) --- .../Polygon_mesh_processing/pmp_compute_normals_test.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp index cd49763c1f7..ef4e8a08f8a 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp @@ -141,10 +141,12 @@ void test(const char* filename) { std::cout << "test " << filename << "..." << std::endl; + // EPECK disabled because it takes too long for the testsuite due to sq roots comparisons, + // but it passed. test_SM(filename); - test_SM(filename); +// test_SM(filename); test_Polyhedron(filename); - test_Polyhedron(filename); +// test_Polyhedron(filename); } int main()