diff --git a/BGL/include/CGAL/boost/graph/named_params_helper.h b/BGL/include/CGAL/boost/graph/named_params_helper.h index 47191997c31..4c432a311c6 100644 --- a/BGL/include/CGAL/boost/graph/named_params_helper.h +++ b/BGL/include/CGAL/boost/graph/named_params_helper.h @@ -112,7 +112,8 @@ namespace CGAL { typedef typename boost::property_traits::value_type type; }; - template + template > class GetVertexPointMap { typedef typename property_map_selector::const_type 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; } 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 10e022ad608..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 @@ -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). @@ -9,88 +9,101 @@ // // // 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 -namespace CGAL{ +#include -namespace Polygon_mesh_processing{ +#include +#include +#include +#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 { 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)); - //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 +void normalize(typename GT::Vector_3& v, const GT& traits) +{ + typedef typename GT::FT FT; - template - typename GT::Vector_3 - triangle_normal(const Point& p0, const Point& p1, const Point& p2 - , const GT& traits) + //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. + const FT norm = CGAL::approximate_sqrt(traits.compute_squared_length_3_object()(v)); + if(norm != FT(0)) { - 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) +{ + 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)); + + //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, FT(1)/FT(2)); +} + +template void sum_normals(const PM& pmesh, typename boost::graph_traits::face_descriptor f, VertexPointMap vpmap, 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 typename boost::property_traits::reference Point_ref; 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)) - { - const Point& pvn = get(vpmap, target(he, pmesh)); - const Point& pvnn = get(vpmap, target(next(he, pmesh), pmesh)); + const Point_ref pv = get(vpmap, v); - Vector n = internal::triangle_normal(pv, pvn, pvnn, traits); + 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)); + + const Vector n = internal::triangle_normal(pv, pvn, pvnn, traits); sum = traits.construct_sum_of_vectors_3_object()(sum, n); he = next(he, pmesh); } } +} // namespace internal /** * \ingroup PMP_normal_grp @@ -124,42 +137,50 @@ 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 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); - sum_normals(pmesh, f, vpmap, normal, traits); + 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)) + if(!traits.equal_3_object()(normal, CGAL::NULL_VECTOR)) internal::normalize(normal, traits); return normal; } +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::parameters::all_default()); +} + /** * \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 @@ -173,22 +194,366 @@ 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, + Face_normal_map face_normals, + 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); + 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 } } +template +void compute_face_normals(const PolygonMesh& pmesh, Face_normal_map face_normals) +{ + 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) +{ + 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 +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, + const FaceNormalVector& face_normals, + const K& traits) +{ + typedef typename K::FT FT; + 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::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(); + for(std::size_t 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(traits.equal_3_object()(ni, nj)) + return ni; + + return traits.construct_sum_of_vectors_3_object()(ni, nj); // not normalized +} + +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; + + 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)) + { + 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); + 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 = cv_3(CGAL::ORIGIN, c); // note that this isn't normalized + } + + return nb; +} + +template +typename GT::Vector_3 +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_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 = cv_3(CGAL::NULL_VECTOR); + + const std::size_t nif = incident_faces.size(); + for(std::size_t i=0; i= 0); + if(sp_bi <= min_sp) + continue; + + if(!does_enclose_other_normals(i, j, -1 /*NA*/, nb, sp_bi, incident_faces, face_normals, traits)) + continue; + + min_sp = sp_bi; + n = nb; + } + } + + return n; +} + +template +typename GT::Vector_3 +compute_most_visible_normal_3_points(const 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; + +#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 = traits.construct_vector_3_object()(CGAL::NULL_VECTOR); + + const std::size_t nif = incident_faces.size(); + for(std::size_t i=0; i(i, j, k, nb, sp_bi, incident_faces, face_normals, traits)) + continue; + + min_sp = sp_bi; + n = nb; + } + } + } + +#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP + std::cout << "Best normal from 3-normals-approach: " << n << std::endl; +#endif + + return n; +} + +// Inspired by Aubry et al. On the most 'normal' normal +template +typename GT::Vector_3 +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::Vector_3 Vector_3; + + std::vector incident_faces; + for(face_descriptor f : CGAL::faces_around_target(halfedge(v, pmesh), pmesh)) + { + if(f == boost::graph_traits::null_face()) + continue; + + incident_faces.push_back(f); + } + + if(incident_faces.size() == 1) + return get(face_normals, incident_faces.front()); + + Vector_3 res = compute_most_visible_normal_2_points(incident_faces, face_normals, traits); + + if(res != CGAL::NULL_VECTOR) // found a valid normal through 2 point min circle + return res; + + CGAL_assertion(incident_faces.size() > 2); + + return compute_most_visible_normal_3_points(incident_faces, face_normals, traits); +} + +template +typename GT::Vector_3 +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::halfedge_descriptor halfedge_descriptor; + + typedef typename GT::FT FT; + typedef typename GT::Vector_3 Vector_3; + typedef typename boost::property_traits::reference Vector_ref; + + 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(); + + 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(!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))); + + //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; +} + +} // end namespace internal + /** * \ingroup PMP_normal_grp * computes the unit normal at vertex `v` as the average of the normals of incident faces. @@ -222,61 +587,101 @@ 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 parameters::choose_parameter; + using parameters::is_default_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 boost::graph_traits::face_descriptor face_descriptor; + + typedef typename GetGeomTraits::type GT; + 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 internal_np::Lookup_named_param_def < - internal_np::face_normal_t, - NamedParameters, - DefaultMap> ::type FaceNormalMap; - FaceNormalMap fnmap = choose_parameter(get_parameter(np, internal_np::face_normal), DefaultMap()); - bool fnmap_valid - = !boost::is_same::value; + 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 boost::graph_traits::halfedge_descriptor halfedge_descriptor; + typedef std::map Face_vector_map; + typedef boost::associative_property_map Default_map; - Vector normal = traits.construct_vector_3_object()(CGAL::NULL_VECTOR); + 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(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; + std::cout << "----------------------------------------------------------------------" << std::endl; + std::cout << "compute vertex at " << get(vpmap, v) + << ", must compute face normals? " << must_compute_face_normals << std::endl; +#endif - halfedge_descriptor he = halfedge(v, pmesh); // handle isolated vertices - if (he==boost::graph_traits::null_halfedge()) return normal; - halfedge_descriptor end = he; - do - { - if (!is_border(he, 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); - } - he = opposite(next(he, pmesh), pmesh); - } while (he != end); + halfedge_descriptor he = halfedge(v, pmesh); + if(he == boost::graph_traits::null_halfedge()) + return CGAL::NULL_VECTOR; - if ( ! traits.equal_3_object()(normal, CGAL::NULL_VECTOR)) + if(must_compute_face_normals) + { + for(face_descriptor f : CGAL::faces_around_target(halfedge(v, pmesh), pmesh)) + { + if(f == boost::graph_traits::null_face()) + continue; + + put(face_normals, f, compute_face_normal(f, pmesh, np)); + } + } + +#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)) << 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 + { +#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); + } + + if(!traits.equal_3_object()(normal, CGAL::NULL_VECTOR)) internal::normalize(normal, traits); + return normal; } +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::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 +* @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 @@ -290,38 +695,82 @@ 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 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)){ - typename Kernel::Vector_3 vec = compute_vertex_normal(v, pmesh, np); - put(vnm, v, vec); + 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::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_map face_normals = choose_parameter(get_parameter(np, internal_np::face_normal), + get(Face_normal_tag(), pmesh)); + const bool must_compute_face_normals = is_default_parameter(get_parameter(np, internal_np::face_normal)); + + if(must_compute_face_normals) + compute_face_normals(pmesh, face_normals, np); + +#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()) + + 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)); + put(vertex_normals, v, n); + +#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 } } +template +void compute_vertex_normals(const PolygonMesh& pmesh, VertexNormalMap vertex_normals) +{ + compute_vertex_normals(pmesh, vertex_normals, CGAL::parameters::all_default()); +} + /** * \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 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 @@ -335,79 +784,27 @@ compute_vertex_normals(const PolygonMesh& pmesh * 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 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)); } -///\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 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()); } -/// \endcond - -} - -} // end of namespace CGAL::Polygon_mesh_processing +} // namespace Polygon_mesh_processing +} // namespace CGAL #endif // CGAL_POLYGON_MESH_PROCESSING_COMPUTE_NORMAL_H 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 ) 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..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 @@ -1,62 +1,165 @@ -#include -#include -#include -#include +// #define CGAL_PMP_COMPUTE_NORMAL_DEBUG -#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_with_sqrt 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 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()); + + 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)); + std::cout.precision(17); + 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; + + std::cout << "Test with Surface_mesh, and kernel: " << typeid(K).name() << std::endl; + + 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; + + std::cout << "Test with Polyhedron, and kernel: " << typeid(K).name() << std::endl; + + 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; + + // 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_Polyhedron(filename); +// test_Polyhedron(filename); } int main() { - test_normals("data/elephant.off"); - test_normals("data/elephant.off"); + std::cout.precision(17); + + 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/U.off"); std::cerr << "All done." << std::endl; return EXIT_SUCCESS;