Merge pull request #4065 from MaelRL/PMP-New_vertex_normal_computations-GF

PMP: add a new way to compute vertex normals

BGL/include/CGAL/boost/graph/named_params_helper.h

@@ -112,7 +112,8 @@ namespace CGAL {
     typedef typename boost::property_traits<PMap>::value_type type;
   };
 
-  template<typename PolygonMesh, typename NamedParameters>
+  template<typename PolygonMesh,
+           typename NamedParameters = Named_function_parameters<bool, internal_np::all_default_t> >
   class GetVertexPointMap
   {
     typedef typename property_map_selector<PolygonMesh, boost::vertex_point_t>::const_type

Polygon_mesh_processing/examples/Polygon_mesh_processing/compute_normals_example.cpp

@@ -6,7 +6,10 @@
 #include <iostream>
 #include <fstream>
 
+namespace PMP = CGAL::Polygon_mesh_processing;
+
 typedef CGAL::Exact_predicates_inexact_constructions_kernel       K;
+
 typedef K::Point_3                                                Point;
 typedef K::Vector_3                                               Vector;
 
@@ -25,24 +28,18 @@ int main(int argc, char* argv[])
     return 1;
   }
 
-  auto fnormals = mesh.add_property_map<face_descriptor, Vector>
+  auto vnormals = mesh.add_property_map<vertex_descriptor, Vector>("v:normals", CGAL::NULL_VECTOR).first;
-      ("f:normals", CGAL::NULL_VECTOR).first;
+  auto fnormals = mesh.add_property_map<face_descriptor, Vector>("f:normals", CGAL::NULL_VECTOR).first;
-  auto vnormals = mesh.add_property_map<vertex_descriptor, Vector>
-      ("v:normals", CGAL::NULL_VECTOR).first;
 
-  CGAL::Polygon_mesh_processing::compute_normals(mesh,
+  PMP::compute_normals(mesh, vnormals, fnormals);
-        vnormals,
+
-        fnormals,
+  std::cout << "Vertex normals :" << std::endl;
-        CGAL::Polygon_mesh_processing::parameters::vertex_point_map(mesh.points()).
+  for(vertex_descriptor vd: vertices(mesh))
-        geom_traits(K()));
+    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;
 }

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,51 +9,62 @@
 // 
 //
 // Author(s)     : Jane Tournois
-
+//                 Mael Rouxel-Labbé
 
 #ifndef CGAL_POLYGON_MESH_PROCESSING_COMPUTE_NORMAL_H
 #define CGAL_POLYGON_MESH_PROCESSING_COMPUTE_NORMAL_H
 
 #include <CGAL/license/Polygon_mesh_processing/Compute_normal.h>
 
-
-#include <CGAL/boost/graph/helpers.h>
-#include <CGAL/boost/graph/properties.h>
-#include <boost/graph/graph_traits.hpp>
-#include <CGAL/Origin.h>
-#include <CGAL/Kernel/global_functions_3.h>
-#include <CGAL/Kernel_traits.h>
-
 #include <CGAL/Polygon_mesh_processing/internal/named_function_params.h>
 #include <CGAL/Polygon_mesh_processing/internal/named_params_helper.h>
 
-#include <boost/type_traits.hpp>
+#include <CGAL/boost/graph/helpers.h>
+#include <CGAL/boost/graph/properties.h>
+#include <CGAL/Dynamic_property_map.h>
+#include <CGAL/Origin.h>
+
+#include <boost/graph/graph_traits.hpp>
+
+#include <iostream>
+#include <utility>
+#include <vector>
+
+#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 <fstream>
+#include <CGAL/Polygon_mesh_processing/bbox.h>
+#endif
 
 namespace CGAL {
-
 namespace Polygon_mesh_processing {
-
 namespace internal {
 
 template <class GT>
 void normalize(typename GT::Vector_3& v, const GT& traits)
 {
-    typename GT::FT norm = CGAL::approximate_sqrt(
+  typedef typename GT::FT                                             FT;
-        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)
+  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);
   }
 }
 
-  template<typename Point
+template<typename Point, typename GT>
-         , typename GT>
 typename GT::Vector_3
-  triangle_normal(const Point& p0, const Point& p1, const Point& p2
+triangle_normal(const Point& p0, const Point& p1, const Point& p2, const GT& traits)
-                , 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));
@@ -61,12 +72,10 @@ namespace internal {
   //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<typename Point, typename PM, typename VertexPointMap, typename Vector
+template<typename Point, typename PM, typename VertexPointMap, typename Vector, typename GT>
-       , typename GT>
 void sum_normals(const PM& pmesh,
                  typename boost::graph_traits<PM>::face_descriptor f,
                  VertexPointMap vpmap,
@@ -76,21 +85,25 @@ void sum_normals(const PM& pmesh,
   typedef typename boost::graph_traits<PM>::vertex_descriptor           vertex_descriptor;
   typedef typename boost::graph_traits<PM>::halfedge_descriptor         halfedge_descriptor;
 
+  typedef typename boost::property_traits<VertexPointMap>::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_ref pvn = get(vpmap, target(he, pmesh));
-    const Point& pvnn = get(vpmap, target(next(he, pmesh), 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);
   }
 }
 
+} // namespace internal
 
 /**
 * \ingroup PMP_normal_grp
@@ -124,9 +137,9 @@ Vector_3
 #else
 typename GetGeomTraits<PolygonMesh, NamedParameters>::type::Vector_3
 #endif
-compute_face_normal(typename boost::graph_traits<PolygonMesh>::face_descriptor f
+compute_face_normal(typename boost::graph_traits<PolygonMesh>::face_descriptor f,
-                    , const PolygonMesh& pmesh
+                    const PolygonMesh& pmesh,
-                    , const NamedParameters& np)
+                    const NamedParameters& np)
 {
   using parameters::choose_parameter;
   using parameters::get_parameter;
@@ -139,10 +152,10 @@ compute_face_normal(typename boost::graph_traits<PolygonMesh>::face_descriptor f
                              get_const_property_map(vertex_point, pmesh));
 
   typedef typename GT::Point_3                                                     Point;
-  typedef typename GT::Vector_3 Vector;
+  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);
-  sum_normals<Point>(pmesh, f, vpmap, normal, traits);
+  internal::sum_normals<Point>(pmesh, f, vpmap, normal, traits);
 
   if(!traits.equal_3_object()(normal, CGAL::NULL_VECTOR))
     internal::normalize(normal, traits);
@@ -150,16 +163,24 @@ compute_face_normal(typename boost::graph_traits<PolygonMesh>::face_descriptor f
   return normal;
 }
 
+template <typename PolygonMesh>
+typename GetGeomTraits<PolygonMesh>::type::Vector_3
+compute_face_normal(typename boost::graph_traits<PolygonMesh>::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<PolygonMesh>::%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<PolygonMesh>::face_descriptor f
 * If `Kernel::FT` does not have a `sqrt()` operation, the square root computation
 * will be done approximately.
 */
-template <typename PolygonMesh
+template <typename PolygonMesh, typename Face_normal_map, typename NamedParameters>
-          , typename FaceNormalMap
+void compute_face_normals(const PolygonMesh& pmesh,
-          , typename NamedParameters>
+                          Face_normal_map face_normals,
-void
+                          const NamedParameters& np)
-compute_face_normals(const PolygonMesh& pmesh
-                   , FaceNormalMap fnm
-                   , const NamedParameters& np)
 {
   typedef typename GetGeomTraits<PolygonMesh,NamedParameters>::type Kernel;
 
-  for(typename boost::graph_traits<PolygonMesh>::face_descriptor f : faces(pmesh)){
+  for(typename boost::graph_traits<PolygonMesh>::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 <typename PolygonMesh, typename Face_normal_map>
+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 <typename GT>
+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 <typename PolygonMesh, typename FaceNormalVector, typename 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<typename boost::graph_traits<PolygonMesh>::face_descriptor>& incident_faces,
+                                const FaceNormalVector& face_normals,
+                                const K& traits)
+{
+  typedef typename K::FT                                                   FT;
+  typedef typename boost::property_traits<FaceNormalVector>::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<nif; ++l)
+  {
+    if(l == i || l == j || l == k)
+      continue;
+
+    const Vector_ref nl = get(face_normals, incident_faces[l]);
+
+    // this is a bound on how much the scalar product between (v1,v2) and (v1, v3) can change
+    // when the angle changes theta_bound := 0.01°
+    // The weird number is thus := max_(theta_i, theta_j)[abs(std::cos(theta_i), std::cos(theta_j))]
+    // with theta_j - theta_i = theta_bound
+    const FT sp_diff_bound = nbn * 0.00017453292431333;
+    const FT sp_bl = sp(nb, nl);
+    if(CGAL::abs(sp_bi - sp_bl) <= sp_diff_bound)
+      continue;
+
+    if(sp_bl < sp_bi)
+      return false;
+  }
+
+  return true;
+}
+
+template <typename GT>
+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>
+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 PolygonMesh, typename FaceNormalVector, typename GT>
+typename GT::Vector_3
+compute_most_visible_normal_2_points(std::vector<typename boost::graph_traits<PolygonMesh>::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<FaceNormalVector>::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<nif; ++i)
+  {
+    for(std::size_t j=i+1; j<nif; ++j)
+    {
+      const Vector_ref ni = get(face_normals, incident_faces[i]);
+      const Vector_ref nj = get(face_normals, incident_faces[j]);
+
+      const Vector_3 nb = compute_normals_bisector(ni, nj, traits);
+
+      // Degeneracies like ni == -nj or a numerical error in the construction of 'nb' can happen.
+      if(traits.equal_3_object()(nb, CGAL::NULL_VECTOR))
+        return CGAL::NULL_VECTOR;
+
+      const FT sp_bi = sp_3(nb, ni);
+
+      CGAL_assertion(sp_bi >= 0);
+      if(sp_bi <= min_sp)
+        continue;
+
+      if(!does_enclose_other_normals<PolygonMesh>(i, j, -1 /*NA*/, nb, sp_bi, incident_faces, face_normals, traits))
+        continue;
+
+      min_sp = sp_bi;
+      n = nb;
+    }
+  }
+
+  return n;
+}
+
+template <typename PolygonMesh, typename FaceNormalVector, typename GT>
+typename GT::Vector_3
+compute_most_visible_normal_3_points(const std::vector<typename boost::graph_traits<PolygonMesh>::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<FaceNormalVector>::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<nif; ++i)
+  {
+    for(std::size_t j=i+1; j<nif; ++j)
+    {
+      for(std::size_t k=j+1; k<nif; ++k)
+      {
+        const Vector_ref ni = get(face_normals, incident_faces[i]);
+        const Vector_ref nj = get(face_normals, incident_faces[j]);
+        const Vector_ref nk = get(face_normals, incident_faces[k]);
+
+        Vector_3 nb = compute_normals_bisector(ni, nj, nk, traits);
+        if(traits.equal_3_object()(nb, CGAL::NULL_VECTOR))
+          return nb;
+
+        FT sp_bi = traits.compute_scalar_product_3_object()(nb, ni);
+        if(sp_bi < FT(0))
+        {
+          nb = traits.construct_opposite_vector_3_object()(nb);
+          sp_bi = - sp_bi;
+        }
+
+        if(sp_bi <= min_sp)
+          continue;
+
+        if(!does_enclose_other_normals<PolygonMesh>(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 PolygonMesh, typename FaceNormalVector, typename GT>
+typename GT::Vector_3
+compute_vertex_normal_most_visible_min_circle(typename boost::graph_traits<PolygonMesh>::vertex_descriptor v,
+                                              const FaceNormalVector& face_normals,
+                                              const PolygonMesh& pmesh,
+                                              const GT& traits)
+{
+  typedef typename boost::graph_traits<PolygonMesh>::face_descriptor       face_descriptor;
+
+  typedef typename GT::Vector_3                                            Vector_3;
+
+  std::vector<face_descriptor> incident_faces;
+  for(face_descriptor f : CGAL::faces_around_target(halfedge(v, pmesh), pmesh))
+  {
+    if(f == boost::graph_traits<PolygonMesh>::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<PolygonMesh>(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<PolygonMesh>(incident_faces, face_normals, traits);
+}
+
+template <typename PolygonMesh, typename FaceNormalVector, typename VertexPointMap, typename GT>
+typename GT::Vector_3
+compute_vertex_normal_as_sum_of_weighted_normals(typename boost::graph_traits<PolygonMesh>::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<PolygonMesh>::halfedge_descriptor      halfedge_descriptor;
+
+  typedef typename GT::FT                                                     FT;
+  typedef typename GT::Vector_3                                               Vector_3;
+  typedef typename boost::property_traits<FaceNormalVector>::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<PolygonMesh>::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<PolygonMesh, NamedParameters>::type::Vector_3
 #endif
 compute_vertex_normal(typename boost::graph_traits<PolygonMesh>::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 boost::graph_traits<PolygonMesh>::halfedge_descriptor      halfedge_descriptor;
+  typedef typename boost::graph_traits<PolygonMesh>::face_descriptor          face_descriptor;
+
   typedef typename GetGeomTraits<PolygonMesh, NamedParameters>::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<PolygonMesh, NamedParameters>::NoMap DefaultMap;
+  typedef typename GetVertexPointMap<PolygonMesh, NamedParameters>::const_type VPMap;
-  typedef typename internal_np::Lookup_named_param_def <
+  VPMap vpmap = choose_parameter(get_parameter(np, internal_np::vertex_point),
-    internal_np::face_normal_t,
+                             get_const_property_map(vertex_point, pmesh));
+
+  typedef std::map<face_descriptor, Vector_3>                                 Face_vector_map;
+  typedef boost::associative_property_map<Face_vector_map>                    Default_map;
+
+  typedef typename internal_np::Lookup_named_param_def<internal_np::face_normal_t,
                                                        NamedParameters,
-    DefaultMap> ::type FaceNormalMap;
+                                                       Default_map>::type     Face_normal_map;
-  FaceNormalMap fnmap = choose_parameter(get_parameter(np, internal_np::face_normal), DefaultMap());
+  Face_vector_map default_fvmap;
-  bool fnmap_valid
+  Face_normal_map face_normals = choose_parameter(get_parameter(np, internal_np::face_normal),
-    = !boost::is_same<FaceNormalMap,
+                                                  Default_map(default_fvmap));
-                      DefaultMap
+  const bool must_compute_face_normals = is_default_parameter(get_parameter(np, internal_np::face_normal));
-                     >::value;
 
-  typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
+#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
 
-  Vector normal = traits.construct_vector_3_object()(CGAL::NULL_VECTOR);
-
-  halfedge_descriptor he = halfedge(v, pmesh);
   // handle isolated vertices
-  if (he==boost::graph_traits<PolygonMesh>::null_halfedge()) return normal;
+  halfedge_descriptor he = halfedge(v, pmesh);
-  halfedge_descriptor end = he;
+  if(he == boost::graph_traits<PolygonMesh>::null_halfedge())
-  do
+    return CGAL::NULL_VECTOR;
+
+  if(must_compute_face_normals)
   {
-    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))
+      if(f == boost::graph_traits<PolygonMesh>::null_face())
-                             : compute_face_normal(face(he, pmesh), pmesh, np);
+        continue;
-      normal = traits.construct_sum_of_vectors_3_object()(normal, n);
+
+      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);
   }
-    he = opposite(next(he, pmesh), pmesh);
-  } while (he != end);
 
   if(!traits.equal_3_object()(normal, CGAL::NULL_VECTOR))
     internal::normalize(normal, traits);
+
   return normal;
 }
 
+template <typename PolygonMesh>
+typename GetGeomTraits<PolygonMesh>::type::Vector_3
+compute_vertex_normal(typename boost::graph_traits<PolygonMesh>::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<PolygonMesh>::%vertex_descriptor` as key type and
+*                         `boost::graph_traits<PolygonMesh>::%vertex_descriptor` as key type and
-    the return type of `compute_vertex_normal()` as value type.
+*                         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<PolygonMesh>::vertex_descript
 * If `Kernel::FT` does not have a `sqrt()` operation, the square root computation
 * will be done approximately.
 */
-template <typename PolygonMesh
+template <typename PolygonMesh, typename VertexNormalMap, typename NamedParameters>
-          , typename VertexNormalMap
+void compute_vertex_normals(const PolygonMesh& pmesh,
-          , typename NamedParameters
+                            VertexNormalMap vertex_normals,
-          >
+                            const NamedParameters& np)
-void
-compute_vertex_normals(const PolygonMesh& pmesh
-                      , VertexNormalMap vnm
-                      , const NamedParameters& np
-                      )
 {
-  typedef typename GetGeomTraits<PolygonMesh,NamedParameters>::type Kernel;
+  using parameters::choose_parameter;
+  using parameters::is_default_parameter;
+  using parameters::get_parameter;
 
-  for(typename boost::graph_traits<PolygonMesh>::vertex_descriptor v : vertices(pmesh)){
+  typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor           vertex_descriptor;
-    typename Kernel::Vector_3 vec = compute_vertex_normal(v, pmesh, np);
+
-    put(vnm, v, vec);
+  typedef typename GetGeomTraits<PolygonMesh,NamedParameters>::type              GT;
+  typedef typename GT::Vector_3                                                  Vector_3;
+
+  typedef CGAL::dynamic_face_property_t<Vector_3>                                Face_normal_tag;
+  typedef typename boost::property_map<PolygonMesh, Face_normal_tag>::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<PolygonMesh, NamedParameters>::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<internal_np::face_normal_t,
+                                                       NamedParameters,
+                                                       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 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 <typename PolygonMesh, typename VertexNormalMap>
+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<PolygonMesh>::%vertex_descriptor` as key type and
+*    `boost::graph_traits<PolygonMesh>::%vertex_descriptor` as key type and
-    `Kernel::Vector_3` as value type.
+*    `Kernel::Vector_3` as value type.
 * @tparam FaceNormalMap a model of `ReadWritePropertyMap` with
-    `boost::graph_traits<PolygonMesh>::%face_descriptor` as key type and
+*    `boost::graph_traits<PolygonMesh>::%face_descriptor` as key type and
-    `Kernel::Vector_3` as value type.
+*    `Kernel::Vector_3` as value type.
 *
 * @param pmesh the polygon mesh
-* @param vnm the property map in which the vertex normals are written
+* @param vertex_normals the property map in which the vertex normals are written
-* @param fnm the property map in which the face 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 <typename PolygonMesh
+template <typename PolygonMesh,
-          , typename VertexNormalMap
+          typename VertexNormalMap, typename FaceNormalMap,
-          , typename FaceNormalMap
+          typename NamedParameters>
-          , typename NamedParameters
+void compute_normals(const PolygonMesh& pmesh,
-          >
+                     VertexNormalMap vertex_normals,
-void
+                     FaceNormalMap face_normals,
-compute_normals(const PolygonMesh& pmesh
+                     const NamedParameters& np)
-                , VertexNormalMap vnm
-                , FaceNormalMap fnm
-                , const NamedParameters& np
-                )
 {
-  compute_face_normals(pmesh, fnm, np);
+  compute_face_normals(pmesh, face_normals, np);
-  compute_vertex_normals(pmesh, vnm, np.face_normal_map(fnm));
+  compute_vertex_normals(pmesh, vertex_normals, np.face_normal_map(face_normals));
 }
 
-///\cond SKIP_IN_MANUAL
-// compute_vertex_normal overloads
-template <typename PolygonMesh>
-typename CGAL::Kernel_traits< typename property_map_value<PolygonMesh, CGAL::vertex_point_t>::type>::Kernel::Vector_3
-compute_vertex_normal(
-  typename boost::graph_traits<PolygonMesh>::vertex_descriptor v,
-  const PolygonMesh& pmesh)
-{
-  return compute_vertex_normal(v, pmesh,
-    CGAL::Polygon_mesh_processing::parameters::all_default());
-}
-
-// compute_vertex_normals overloads
-template <typename PolygonMesh, typename VertexNormalMap>
-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 PolygonMesh>
-typename CGAL::Kernel_traits < typename property_map_value<PolygonMesh, CGAL::vertex_point_t>::type>::Kernel::Vector_3
-compute_face_normal(
-  typename boost::graph_traits<PolygonMesh>::face_descriptor f,
-  const PolygonMesh& pmesh)
-{
-  return compute_face_normal(f, pmesh,
-    CGAL::Polygon_mesh_processing::parameters::all_default());
-}
-
-// compute_face_normals overload
-template <typename PolygonMesh, typename FaceNormalMap>
-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 <typename PolygonMesh, typename VertexNormalMap, typename FaceNormalMap>
-void
+void compute_normals(const PolygonMesh& pmesh,
-compute_normals(const PolygonMesh& pmesh,
+                     VertexNormalMap vertex_normals,
-                VertexNormalMap vnm,
+                     FaceNormalMap face_normals)
-                FaceNormalMap fnm)
 {
-  compute_normals(pmesh, vnm, fnm,
+  compute_normals(pmesh, vertex_normals, face_normals, CGAL::parameters::all_default());
-    CGAL::Polygon_mesh_processing::parameters::all_default());
 }
 
-/// \endcond
+} // namespace Polygon_mesh_processing
-
+} // namespace CGAL
-}
-
-} // end of namespace CGAL::Polygon_mesh_processing
 
 #endif // CGAL_POLYGON_MESH_PROCESSING_COMPUTE_NORMAL_H

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 )
 

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

Polygon_mesh_processing/test/Polygon_mesh_processing/pmp_compute_normals_test.cpp

@@ -1,62 +1,165 @@
-#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+// #define CGAL_PMP_COMPUTE_NORMAL_DEBUG
-#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
-#include <CGAL/Surface_mesh.h>
-#include <CGAL/Polygon_mesh_processing/compute_normal.h>
 
-#include <CGAL/Timer.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel_with_sqrt.h>
+
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/Polyhedron_3.h>
+
+#include <CGAL/Polygon_mesh_processing/compute_normal.h>
 
 #include <iostream>
 #include <fstream>
 
-typedef CGAL::Exact_predicates_inexact_constructions_kernel Epic;
+typedef CGAL::Exact_predicates_inexact_constructions_kernel          EPICK;
-typedef CGAL::Exact_predicates_exact_constructions_kernel Epec;
+typedef CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt  EPECK;
+
+typedef CGAL::Surface_mesh<EPICK::Point_3>                           EPICK_SM;
+typedef CGAL::Surface_mesh<EPECK::Point_3>                           EPECK_SM;
+
+namespace PMP = CGAL::Polygon_mesh_processing;
+
+template <typename K, typename Mesh, typename VertexNormalPmap, typename FaceNormalPmap>
+void test(const Mesh& mesh,
+          const VertexNormalPmap& vnormals,
+          const FaceNormalPmap& fnormals)
+{
+  typedef typename K::Vector_3                                            Vector;
+
+  typedef typename boost::graph_traits<Mesh>::vertex_descriptor           vertex_descriptor;
+  typedef typename boost::graph_traits<Mesh>::face_descriptor             face_descriptor;
+
+  typedef typename CGAL::GetVertexPointMap<Mesh>::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<typename K>
-void test_normals(const char* file_name)
+void test_SM(const char* file_name)
 {
-  typedef typename K::Point_3 Point;
+  typedef CGAL::Surface_mesh<typename K::Point_3>                         SM;
-  typedef typename K::Vector_3 Vector;
+  typedef typename boost::graph_traits<SM>::vertex_descriptor             vertex_descriptor;
-  typedef CGAL::Surface_mesh<Point> Surface_mesh;
+  typedef typename boost::graph_traits<SM>::face_descriptor               face_descriptor;
-  typedef typename boost::graph_traits<Surface_mesh>::vertex_descriptor vertex_descriptor;
-  typedef typename boost::graph_traits<Surface_mesh>::face_descriptor face_descriptor;
 
-  Surface_mesh mesh;
+  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)){
+  if(!(input >> mesh))
-    std::cerr << "Error: cannot read Surface_mesh : " << file_name << "\n";
+  {
+    std::cerr << "Error: cannot read " << file_name << " as a Surface_mesh\n";
     assert(false);
   }
 
-  typename Surface_mesh::template Property_map<face_descriptor, Vector> fnormals;
+  typename SM::template Property_map<vertex_descriptor, Vector> vnormals;
-  bool created;
+  vnormals = mesh.template add_property_map<vertex_descriptor, Vector>("v:normals", CGAL::NULL_VECTOR).first;
-  boost::tie(fnormals, created) = mesh.template add_property_map<face_descriptor,Vector>("f:normals",Vector(0,0,0));
+  typename SM::template Property_map<face_descriptor, Vector> fnormals;
-  CGAL::Polygon_mesh_processing::compute_face_normals(mesh, fnormals);
+  fnormals = mesh.template add_property_map<face_descriptor, Vector>("f:normals", CGAL::NULL_VECTOR).first;
-  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<vertex_descriptor,Vector> vnormals;
+  test<K>(mesh, vnormals, fnormals);
+}
 
-  boost::tie(vnormals, created) = mesh.template add_property_map<vertex_descriptor,Vector>("v:normals",Vector(0,0,0));
+template<typename K>
-  CGAL::Polygon_mesh_processing::compute_vertex_normals(mesh, vnormals);
+void test_Polyhedron(const char* file_name)
-  CGAL::Polygon_mesh_processing::compute_vertex_normals(mesh, vnormals,
+{
-    CGAL::Polygon_mesh_processing::parameters::vertex_point_map(mesh.points())); 
+  typedef CGAL::Polyhedron_3<K>                                           Polyhedron;
-  CGAL::Polygon_mesh_processing::compute_vertex_normals(mesh, vnormals,
+  typedef typename boost::graph_traits<Polyhedron>::vertex_descriptor     vertex_descriptor;
-    CGAL::Polygon_mesh_processing::parameters::vertex_point_map(mesh.points()).geom_traits(K()));
+  typedef typename boost::graph_traits<Polyhedron>::face_descriptor       face_descriptor;
 
-  CGAL::Polygon_mesh_processing::compute_normals(mesh, vnormals, fnormals);
+  typedef typename K::Vector_3                                            Vector;
-  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()));
 
+  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<vertex_descriptor, Vector>                             VertexNormalMap;
+  typedef boost::associative_property_map<VertexNormalMap>                VertexNormalPMap;
+  typedef std::map<face_descriptor, Vector>                               FaceNormalMap;
+  typedef boost::associative_property_map<FaceNormalMap>                  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<K>(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<EPICK>(filename);
+//  test_SM<EPECK>(filename);
+  test_Polyhedron<EPICK>(filename);
+//  test_Polyhedron<EPECK>(filename);
 }
 
 int main()
 {
-  test_normals<Epic>("data/elephant.off");
+  std::cout.precision(17);
-  test_normals<Epec>("data/elephant.off");
+
+  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;