Clean code

Surface_mesh_simplification/examples/Surface_mesh_simplification/gh_statistics.cpp

@@ -284,7 +284,7 @@ void hausdorff_errors_dir(const fs::path& dir, const fs::path& out, double ratio
 
 void write_aspect_ratios(fs::ofstream& out, const Surface_mesh& mesh) 
 {
-  for (auto face : mesh.faces())
+  for(auto face : faces(mesh))
   {
     out << std::to_string(PMP::face_aspect_ratio(face, mesh)) << '\n'; 
   }

Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_policies.h

@@ -24,10 +24,6 @@
 namespace CGAL {
 namespace Surface_mesh_simplification {
 
-// derived class implements functions used in the base class that
-// takes the derived class as template argument - see "CRTP"
-//
-// implements classic plane quadrics
 template<typename TriangleMesh, typename GeomTraits>
 class GarlandHeckbert_policies :
   public internal::GarlandHeckbert_placement_base<
@@ -44,16 +40,13 @@ class GarlandHeckbert_policies :
              CGAL::dynamic_vertex_property_t<Eigen::Matrix<typename GeomTraits::FT, 4, 4, Eigen::DontAlign> >
            >::type,
            GeomTraits,
-    GarlandHeckbert_policies<TriangleMesh, GeomTraits>
+           GarlandHeckbert_policies<TriangleMesh, GeomTraits> >
-  >
 {
-
 public:
   typedef typename GeomTraits::FT                                              FT;
 
-  // This is ugly, we later only use the Mat_4 from the 
+  // This is ugly, we later only use the Mat_4 from the Cost_base, but we want to define
-  // Cost_base, but we want to define the matrix here already so it's nicer to define
+  // the matrix here already so it's nicer to define Cost_base in the first place
-  // Cost_base in the first place
   typedef typename Eigen::Matrix<FT, 4, 4, Eigen::DontAlign>                   GH_matrix;
   typedef CGAL::dynamic_vertex_property_t<GH_matrix>                           Cost_property;
 
@@ -67,54 +60,57 @@ class GarlandHeckbert_policies :
     Vertex_cost_map, GeomTraits, GarlandHeckbert_policies<TriangleMesh, GeomTraits>
     >                                                                          Cost_base;
 
-  // both types are the same, this is so we avoid casting back to the base class in
+  // both types are the same, this is so we avoid casting back to the base class in get_cost() or get_placement()
-  // get_cost() or get_placement()
   typedef GarlandHeckbert_policies                                             Get_cost;
   typedef GarlandHeckbert_policies                                             Get_placement;
 
-  // introduce both operators into scope so we get an overload operator()
+  // these 'using' directives are needed to choose between the definitions of these types
-  // this is needed since Get_cost and Get_placement are the same type
-  using Cost_base::operator();
-  using Placement_base::operator();
-
-  // these using directives are needed to choose between the definitions of these types
   // in Cost_base and Placement_base (even though they are the same)
   using typename Cost_base::Mat_4;
   using typename Cost_base::Col_4;
   using typename Cost_base::Point_3;
   using typename Cost_base::Vector_3;
 
-  GarlandHeckbert_policies(TriangleMesh& tmesh, FT dm = FT(100)) 
+private:
+  Vertex_cost_map vcm_;
+
+public:
+  GarlandHeckbert_policies(TriangleMesh& tmesh,
+                           FT dm = FT(100)) // unused @tmp
   {
     vcm_ = get(Cost_property(), tmesh);
 
-    /**
+    // initialize the two base class cost matrices (protected members)
-     * initialize the two base class cost matrices (protected members)
-     */
     Cost_base::init_vcm(vcm_);
     Placement_base::init_vcm(vcm_);
   }
 
+  const Get_cost& get_cost() const { return *this; }
+  const Get_placement& get_placement() const { return *this; }
+  Vertex_cost_map get_vcm() const { return vcm_; }
+
+  // introduce both operators into scope so we get an overload operator()
+  // this is needed since Get_cost and Get_placement are the same type
+  using Cost_base::operator();
+  using Placement_base::operator();
+
+public:
   template<typename VertexPointMap>
-  Mat_4 construct_quadric_from_face(
+  Mat_4 construct_quadric_from_face(typename boost::graph_traits<TriangleMesh>::face_descriptor f,
-      const VertexPointMap& point_map,
                                     const TriangleMesh& tmesh,
-      typename boost::graph_traits<TriangleMesh>::face_descriptor f,
+                                    const VertexPointMap point_map,
                                     const GeomTraits& gt) const
   {
-    return internal::construct_classic_plane_quadric_from_face(
+    return internal::construct_classic_plane_quadric_from_face(f, tmesh, point_map, gt);
-        point_map, tmesh, f, gt);
   }
 
   template<typename VertexPointMap>
-  Mat_4 construct_quadric_from_edge(
+  Mat_4 construct_quadric_from_edge(typename boost::graph_traits<TriangleMesh>::halfedge_descriptor he,
-      const VertexPointMap& point_map,
                                     const TriangleMesh& tmesh,
-      typename boost::graph_traits<TriangleMesh>::halfedge_descriptor he,
+                                    const VertexPointMap point_map,
                                     const GeomTraits& gt) const
   {
-    return internal::construct_classic_plane_quadric_from_edge(
+    return internal::construct_classic_plane_quadric_from_edge(he, tmesh, point_map, gt);
-        point_map, tmesh, he, gt);
   }
 
   Col_4 construct_optimal_point(const Mat_4& quadric, const Col_4& p0, const Col_4& p1) const
@@ -122,21 +118,10 @@ class GarlandHeckbert_policies :
     return internal::construct_optimal_point_singular<GeomTraits>(quadric, p0, p1);
   }
 
-  const Get_cost& get_cost() const { return *this; }
+  Mat_4 construct_quadric_from_normal(const Vector_3& normal,
-  const Get_placement& get_placement() const { return *this; }
+                                      const Point_3& point,
-
-  // return a copy of the vertex cost map to inspect quadrics
-  Vertex_cost_map get_vcm() const { return vcm_; }
-  
-  private:
-  Vertex_cost_map vcm_;
-
-  // helper function to construct quadrics
-
-  Mat_4 construct_quadric_from_normal(const Vector_3& normal, const Point_3& point, 
                                       const GeomTraits& gt) const
   {
-
     auto dot_product = gt.compute_scalar_product_3_object();
     auto construct_vector = gt.construct_vector_3_object();
 

Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_probabilistic_policies.h

@@ -14,55 +14,49 @@
 #define CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_GARLANDHECKBERT_PROBABILISTIC_POLICIES_H
 
 #include <CGAL/license/Surface_mesh_simplification.h>
+
 #include <CGAL/Surface_mesh_simplification/internal/Common.h>
 #include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/GarlandHeckbert_policy_base.h>
 #include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/GarlandHeckbert_functions.h>
+
 #include <Eigen/Dense>
 
 namespace CGAL {
 namespace Surface_mesh_simplification {
 
-// derived class implements functions used in the base class that
+// Implements probabilistic plane quadrics,
-// takes the derived class as template argument - see "CRTP"
+// optionally takes a face variance map giving a per-face variance
-//
+template<typename TriangleMesh,
-// implements probabilistic plane quadrics, optionally takes a face variance map
+         typename GeomTraits,
-// giving a per-face variance 
+         typename FaceVarianceMap = Constant_property_map<
-template<typename TriangleMesh, typename GeomTraits, typename 
-  FaceVarianceMap = Constant_property_map<
                                       typename boost::graph_traits<TriangleMesh>::face_descriptor,
-             std::pair<typename GeomTraits::FT, typename GeomTraits::FT>
+                                      std::pair<typename GeomTraits::FT, typename GeomTraits::FT> > >
-  >> 
+class GarlandHeckbert_probabilistic_policies
-class GarlandHeckbert_probabilistic_policies :
+  : public internal::GarlandHeckbert_placement_base<
-  public internal::GarlandHeckbert_placement_base<
              typename boost::property_map<
                TriangleMesh,
                CGAL::dynamic_vertex_property_t<Eigen::Matrix<typename GeomTraits::FT, 4, 4, Eigen::DontAlign> >
            >::type,
            GeomTraits,
-    GarlandHeckbert_probabilistic_policies<TriangleMesh, GeomTraits, FaceVarianceMap>
+           GarlandHeckbert_probabilistic_policies<TriangleMesh, GeomTraits, FaceVarianceMap> >,
-  >,
     public internal::GarlandHeckbert_cost_base<
              typename boost::property_map<
                TriangleMesh,
                CGAL::dynamic_vertex_property_t<Eigen::Matrix<typename GeomTraits::FT, 4, 4, Eigen::DontAlign>>
              >::type,
              GeomTraits,
-    GarlandHeckbert_probabilistic_policies<TriangleMesh, GeomTraits, FaceVarianceMap>
+             GarlandHeckbert_probabilistic_policies<TriangleMesh, GeomTraits, FaceVarianceMap> >
-  >
 {
-  
   typedef typename GeomTraits::FT                                              FT;
 
-  typedef typename boost::graph_traits<TriangleMesh>::face_descriptor face_descriptor; 
   typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor      halfedge_descriptor;
+  typedef typename boost::graph_traits<TriangleMesh>::face_descriptor          face_descriptor;
 
   typedef typename boost::property_traits<FaceVarianceMap>::value_type         Face_variance;
 
 public:
-
   typedef typename Eigen::Matrix<FT, 4, 4, Eigen::DontAlign>                   GH_matrix;
   typedef CGAL::dynamic_vertex_property_t<GH_matrix>                           Cost_property;
-
   typedef typename boost::property_map<TriangleMesh, Cost_property>::type      Vertex_cost_map;
 
   typedef internal::GarlandHeckbert_placement_base<
@@ -78,11 +72,6 @@ class GarlandHeckbert_probabilistic_policies :
   typedef GarlandHeckbert_probabilistic_policies                               Get_cost;
   typedef GarlandHeckbert_probabilistic_policies                               Get_placement;
 
-  // so that operator() gets overloaded, this is needed because now Get_cost and Get_placement
-  // are the same
-  using Cost_base::operator();
-  using Placement_base::operator();
-
   // these using directives are needed to choose between the definitions of these types
   // in Cost_base and Placement_base (even though they are the same)
   using typename Cost_base::Mat_4;
@@ -90,6 +79,19 @@ class GarlandHeckbert_probabilistic_policies :
   using typename Cost_base::Point_3;
   using typename Cost_base::Vector_3;
 
+private:
+  Vertex_cost_map vcm_;
+
+  // magic number determined by some testing
+  static constexpr FT good_default_variance_unit = 0.05;
+
+  // magic number - for most use cases, there is no input variance, so it makes sense to
+  // set the positional variance to a smaller value than the normal variance
+  static constexpr FT position_variance_factor = 0.1;
+
+  FaceVarianceMap face_variance_map;
+
+public:
   // default discontinuity multiplier is 100
   GarlandHeckbert_probabilistic_policies(TriangleMesh& tmesh)
     : GarlandHeckbert_probabilistic_policies(tmesh, 100)
@@ -107,8 +109,9 @@ class GarlandHeckbert_probabilistic_policies :
 
     // try to initialize the face variance map using the estimated variance
     // parameters are constants defined for this class
-    face_variance_map = FaceVarianceMap { 
+    face_variance_map =
-      internal::estimate_variances(tmesh, GeomTraits(), good_default_variance_unit, 
+      FaceVarianceMap { internal::estimate_variances(tmesh, GeomTraits(),
+                                                     good_default_variance_unit,
                                                      position_variance_factor) };
   }
 
@@ -123,66 +126,52 @@ class GarlandHeckbert_probabilistic_policies :
     Placement_base::init_vcm(vcm_);
   }
 
+  const Get_cost& get_cost() const { return *this; }
+  const Get_placement& get_placement() const { return *this; }
+
+  // so that operator() gets overloaded, this is needed because now Get_cost and Get_placement are the same
+  using Cost_base::operator();
+  using Placement_base::operator();
+
+public:
   template<typename VPM>
-  Mat_4 construct_quadric_from_face(
+  Mat_4 construct_quadric_from_face(face_descriptor f,
-      const VPM& point_map,
                                     const TriangleMesh& tmesh,
-      face_descriptor f, 
+                                    const VPM point_map,
                                     const GeomTraits& gt) const
   {
-    const Vector_3 normal = internal::construct_unit_normal_from_face<
+    const Vector_3 normal = internal::construct_unit_normal_from_face(f, tmesh, point_map, gt);
-      VPM, TriangleMesh, GeomTraits>(point_map, tmesh, f, gt);
-    
     const Point_3 p = get(point_map, source(halfedge(f, tmesh), tmesh));
 
     FT n_variance;
     FT p_variance;
-
     std::tie(n_variance, p_variance) = get(face_variance_map, f);
 
     return internal::construct_prob_plane_quadric_from_normal(normal, p, gt, n_variance, p_variance);
   }
 
   template<typename VPM>
-  Mat_4 construct_quadric_from_edge(
+  Mat_4 construct_quadric_from_edge(halfedge_descriptor he,
-      const VPM& point_map,
                                     const TriangleMesh& tmesh,
-      halfedge_descriptor he, 
+                                    const VPM point_map,
                                     const GeomTraits& gt) const
   {
-    const Vector_3 normal = internal::construct_edge_normal(point_map, tmesh, he, gt);
+    const Vector_3 normal = internal::construct_edge_normal(he, tmesh, point_map, gt);
-
     const Point_3 p = get(point_map, source(he, tmesh));
 
     FT n_variance;
     FT p_variance;
-
     std::tie(n_variance, p_variance) = get(face_variance_map, face(he, tmesh));
 
     return internal::construct_prob_plane_quadric_from_normal(normal, p, gt, n_variance, p_variance);
   }
 
-  Col_4 construct_optimal_point(const Mat_4& quadric, const Col_4& p0, 
+  Col_4 construct_optimal_point(const Mat_4& quadric,
+                                const Col_4& p0,
                                 const Col_4& p1) const
   {
     return internal::construct_optimal_point_invertible<GeomTraits>(quadric);
   }
-  
-  const Get_cost& get_cost() const { return *this; }
-  
-  const Get_placement& get_placement() const { return *this; }
-  
-  private:
-  Vertex_cost_map vcm_;
-  
-  // magic number determined by some testing 
-  static constexpr FT good_default_variance_unit = 0.05;
-
-  // magic number - for most use cases, there is no input variance, so it makes sense to
-  // set the positional variance to a smaller value than the normal variance
-  static constexpr FT position_variance_factor = 0.1;
-    
-  FaceVarianceMap face_variance_map;
 };
 
 } // namespace Surface_mesh_simplification

Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_probabilistic_tri_policies.h

@@ -13,8 +13,10 @@
 #define CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_GARLANDHECKBERT_PROBABILISTIC_TRI_POLICIES_H
 
 #include <CGAL/license/Surface_mesh_simplification.h>
+
 #include <CGAL/Surface_mesh_simplification/internal/Common.h>
 #include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/GarlandHeckbert_policy_base.h>
+
 #include <Eigen/Dense>
 
 #include <CGAL/boost/graph/Named_function_parameters.h>
@@ -28,12 +30,13 @@ namespace Surface_mesh_simplification {
 // derives from cost_base and placement_base
 // implements probabilistic triangle faces and optionally takes a face variance map
 // analogously to probabilistic plane quadrics
-template<typename TriangleMesh, typename GeomTraits, typename
+template<typename TriangleMesh,
-  FaceVarianceMap = Constant_property_map<
+         typename GeomTraits,
-  typename boost::graph_traits<TriangleMesh>::face_descriptor, typename GeomTraits::FT
+         typename FaceVarianceMap =
-   >>
+           Constant_property_map<typename boost::graph_traits<TriangleMesh>::face_descriptor,
-class GarlandHeckbert_probabilistic_tri_policies :
+                                 typename GeomTraits::FT> >
-  public internal::GarlandHeckbert_placement_base<
+class GarlandHeckbert_probabilistic_tri_policies
+  : public internal::GarlandHeckbert_placement_base<
              typename boost::property_map<
                TriangleMesh,
                CGAL::dynamic_vertex_property_t<Eigen::Matrix<typename GeomTraits::FT, 4, 4, Eigen::DontAlign>>
@@ -47,10 +50,8 @@ class GarlandHeckbert_probabilistic_tri_policies :
              CGAL::dynamic_vertex_property_t<Eigen::Matrix<typename GeomTraits::FT, 4, 4, Eigen::DontAlign>>
            >::type,
            GeomTraits,
-    GarlandHeckbert_probabilistic_tri_policies<TriangleMesh, GeomTraits, FaceVarianceMap>
+           GarlandHeckbert_probabilistic_tri_policies<TriangleMesh, GeomTraits, FaceVarianceMap> >
-  >
 {
-
   typedef typename boost::property_traits<FaceVarianceMap>::value_type         Face_variance;
 
 public:
@@ -59,7 +60,6 @@ class GarlandHeckbert_probabilistic_tri_policies :
   typedef typename Eigen::Matrix<FT, 3, 3, Eigen::DontAlign>                   Mat_3;
   typedef typename Eigen::Matrix<FT, 4, 4, Eigen::DontAlign>                   GH_matrix;
   typedef CGAL::dynamic_vertex_property_t<GH_matrix>                           Cost_property;
-
   typedef typename boost::property_map<TriangleMesh, Cost_property>::type      Vertex_cost_map;
 
   typedef internal::GarlandHeckbert_placement_base<
@@ -75,11 +75,6 @@ class GarlandHeckbert_probabilistic_tri_policies :
   typedef GarlandHeckbert_probabilistic_tri_policies                           Get_cost;
   typedef GarlandHeckbert_probabilistic_tri_policies                           Get_placement;
 
-  // so that operator() gets overloaded, this is needed because now Get_cost and Get_placement
-  // are the same
-  using Cost_base::operator();
-  using Placement_base::operator();
-
   // these using directives are needed to choose between the definitions of these types
   // in Cost_base and Placement_base (even though they are the same)
   using typename Cost_base::Mat_4;
@@ -87,6 +82,18 @@ class GarlandHeckbert_probabilistic_tri_policies :
   using typename Cost_base::Point_3;
   using typename Cost_base::Vector_3;
 
+private:
+  Vertex_cost_map vcm_;
+  FaceVarianceMap face_variance_map;
+
+  // same meaning as for probabilistic plane quadrics
+  static constexpr FT good_default_variance_unit = 0.05;
+
+  // this is only used when we fall back to probabilistic planes for the discontinuous edges,
+  // the actual triangle quadric calculation only uses the normal variance
+  static constexpr FT position_variance_factor = 1;
+
+public:
   // default discontinuity multiplier is 100
   GarlandHeckbert_probabilistic_tri_policies(TriangleMesh& tmesh)
     : GarlandHeckbert_probabilistic_tri_policies(tmesh, 100)
@@ -107,6 +114,7 @@ class GarlandHeckbert_probabilistic_tri_policies :
 
     std::tie(variance, discard_position) = internal::estimate_variances(tmesh, GeomTraits(),
         good_default_variance_unit, position_variance_factor);
+
     // see probabilistic plane quadrics
     face_variance_map = FaceVarianceMap { variance };
   }
@@ -116,9 +124,8 @@ class GarlandHeckbert_probabilistic_tri_policies :
                                              const FaceVarianceMap* fvm)
     : Cost_base(dm), face_variance_map(fvm)
   {
-
     // we need positive variances so that we always get an invertible matrix
-    CGAL_precondition(sdn > 0.0 && sdp > 0.0);
+//    CGAL_precondition(sdn > 0. && sdp > 0.); // @fixme what was that, check history
 
     // initialize the private variable vcm so it's lifetime is bound to that of the policy's
     vcm_ = get(Cost_property(), tmesh);
@@ -128,31 +135,35 @@ class GarlandHeckbert_probabilistic_tri_policies :
     Placement_base::init_vcm(vcm_);
   }
 
+  const Get_cost& get_cost() const { return *this; }
+  const Get_placement& get_placement() const { return *this; }
+
+  // so that operator() gets overloaded, this is needed because now Get_cost and Get_placement are the same
+  using Cost_base::operator();
+  using Placement_base::operator();
+
+public:
   template<typename VPM, typename TM>
-  Mat_4 construct_quadric_from_face(
+  Mat_4 construct_quadric_from_face(typename boost::graph_traits<TM>::face_descriptor f,
-      const VPM& point_map,
                                     const TM& tmesh,
-      typename boost::graph_traits<TM>::face_descriptor f, 
+                                    const VPM point_map,
                                     const GeomTraits& gt) const
   {
     FT variance = get(face_variance_map, f);
 
-    return internal::construct_prob_triangle_quadric_from_face(
+    return internal::construct_prob_triangle_quadric_from_face(f, variance, tmesh, point_map, gt);
-        point_map, tmesh, f, variance, gt);
   }
 
   // we don't have a sensible way to construct a triangle quadric
   // from an edge, so we fall back to probabilistic plane quadrics here
   template<typename VPM, typename TM>
-  Mat_4 construct_quadric_from_edge(
+  Mat_4 construct_quadric_from_edge(typename boost::graph_traits<TM>::halfedge_descriptor he,
-      const VPM& point_map,
                                     const TM& tmesh,
-      typename boost::graph_traits<TM>::halfedge_descriptor he, 
+                                    const VPM point_map,
                                     const GeomTraits& gt) const
   {
     // same as in probabilistic plane policy
-    const Vector_3 normal = internal::construct_edge_normal(point_map, tmesh, he, gt);
+    const Vector_3 normal = internal::construct_edge_normal(he, tmesh, point_map, gt);
-
     const Point_3 p = get(point_map, source(he, tmesh));
 
     FT variance = get(face_variance_map, face(he, tmesh));
@@ -165,20 +176,6 @@ class GarlandHeckbert_probabilistic_tri_policies :
   {
     return internal::construct_optimal_point_invertible<GeomTraits>(quadric);
   }
-  const Get_cost& get_cost() const { return *this; }
-  const Get_placement& get_placement() const { return *this; }
-    
-  private:
-  Vertex_cost_map vcm_;
-  
-  FaceVarianceMap face_variance_map;
-  
-  // same meaning as for probabilistic plane quadrics
-  static constexpr FT good_default_variance_unit = 0.05;
-
-  // this is only used when we fall back to probabilistic planes for the discontinuous edges,
-  // the actual triangle quadric calculation only uses the normal variance
-  static constexpr FT position_variance_factor = 1;
 };
 
 } // namespace Surface_mesh_simplification

Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_triangle_policies.h

@@ -18,7 +18,9 @@
 #include <CGAL/Surface_mesh_simplification/internal/Common.h>
 #include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/GarlandHeckbert_policy_base.h>
 #include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/GarlandHeckbert_functions.h>
+
 #include <Eigen/Dense>
+
 #include <iostream>
 
 namespace CGAL {
@@ -29,8 +31,8 @@ namespace Surface_mesh_simplification {
 //
 // implements classic triangle policies
 template<typename TriangleMesh, typename GeomTraits>
-class GarlandHeckbert_triangle_policies : 
+class GarlandHeckbert_triangle_policies
-  public internal::GarlandHeckbert_cost_base<
+  : public internal::GarlandHeckbert_cost_base<
              typename boost::property_map<
                TriangleMesh,
                CGAL::dynamic_vertex_property_t<Eigen::Matrix<typename GeomTraits::FT, 4, 4, Eigen::DontAlign>>
@@ -44,8 +46,7 @@ class GarlandHeckbert_triangle_policies :
              CGAL::dynamic_vertex_property_t<Eigen::Matrix<typename GeomTraits::FT, 4, 4, Eigen::DontAlign>>
            >::type,
            GeomTraits,
-    GarlandHeckbert_triangle_policies<TriangleMesh, GeomTraits>
+           GarlandHeckbert_triangle_policies<TriangleMesh, GeomTraits> >
-  >
 {
 public:
   typedef typename GeomTraits::FT                                              FT;
@@ -64,9 +65,13 @@ class GarlandHeckbert_triangle_policies :
     Vertex_cost_map, GeomTraits, GarlandHeckbert_triangle_policies<TriangleMesh, GeomTraits>
     >                                                                          Cost_base;
 
-  using Cost_base::operator();
+  typedef GarlandHeckbert_triangle_policies                                    Get_cost;
-  using Placement_base::operator();
+  typedef GarlandHeckbert_triangle_policies                                    Get_placement;
 
+private:
+  Vertex_cost_map vcm_;
+
+public:
   GarlandHeckbert_triangle_policies(TriangleMesh& tmesh, FT dm = FT(100))
   {
     vcm_ = get(Cost_property(), tmesh);
@@ -75,42 +80,36 @@ class GarlandHeckbert_triangle_policies :
     Placement_base::init_vcm(vcm_);
   }
 
+public:
+  const Get_cost& get_cost() const { return *this; }
+  const Get_placement& get_placement() const { return *this; }
+
+  using Cost_base::operator();
+  using Placement_base::operator();
+
+public:
   template<typename VertexPointMap>
-  Mat_4 construct_quadric_from_face(
+  Mat_4 construct_quadric_from_face(typename boost::graph_traits<TriangleMesh>::face_descriptor f,
-      const VertexPointMap& point_map,
                                     const TriangleMesh& tmesh,
-      typename boost::graph_traits<TriangleMesh>::face_descriptor f,
+                                    const VertexPointMap point_map,
                                     const GeomTraits& gt) const
   {
-    return internal::construct_classic_triangle_quadric_from_face(
+    return internal::construct_classic_triangle_quadric_from_face(f, tmesh, point_map, gt);
-        point_map, tmesh, f, gt);
   }
 
   template<typename VertexPointMap>
-  Mat_4 construct_quadric_from_edge(
+  Mat_4 construct_quadric_from_edge(typename boost::graph_traits<TriangleMesh>::halfedge_descriptor he,
-      const VertexPointMap& point_map,
                                     const TriangleMesh& tmesh,
-      typename boost::graph_traits<TriangleMesh>::halfedge_descriptor he,
+                                    const VertexPointMap point_map,
                                     const GeomTraits& gt) const
   {
-    return internal::construct_classic_plane_quadric_from_edge(
+    return internal::construct_classic_plane_quadric_from_edge(he, tmesh, point_map, gt);
-        point_map, tmesh, he, gt);
   }
 
   Col_4 construct_optimal_point(const Mat_4& quadric, const Col_4& p0, const Col_4& p1) const
   {
     return internal::construct_optimal_point_singular<GeomTraits>(quadric, p0, p1);
   }
-  
-  typedef GarlandHeckbert_triangle_policies Get_cost;
-  typedef GarlandHeckbert_triangle_policies Get_placement;
-  
-  const Get_cost& get_cost() const { return *this; }
-  const Get_placement& get_placement() const { return *this; }
-  
-  private:
-  
-  Vertex_cost_map vcm_;
 };
 
 } // namespace Surface_mesh_simplification

Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/GarlandHeckbert_functions.h

@@ -17,6 +17,7 @@
 #include <CGAL/Surface_mesh_simplification/internal/Common.h>
 
 #include <Eigen/Dense>
+
 #include <iostream>
 
 namespace CGAL {
@@ -54,10 +55,8 @@ bool invert_matrix_4(const Matrix& m, Matrix& im)
       - m(0, 3) * ( m(1, 0) * A1223 - m(1, 1) * A0223 + m(1, 2) * A0123 );
   det = 1 / det;
 
-  if (det == 0.0)
+  if(det == 0.)
-  {
     return false;
-  }
 
   // we never actually use values other than those in the third column,
   // so might as well not calculate them
@@ -82,8 +81,8 @@ bool invert_matrix_4(const Matrix& m, Matrix& im)
 }
 
 template<typename GeomTraits>
-Eigen::Matrix<typename GeomTraits::FT, 3, 1> vector_to_col_vec(
+Eigen::Matrix<typename GeomTraits::FT, 3, 1>
-    const typename GeomTraits::Vector_3& v) 
+vector_to_col_vec(const typename GeomTraits::Vector_3& v)
 {
   Eigen::Matrix<typename GeomTraits::FT, 3, 1> col {v.x(), v.y(), v.z()};
   return col;
@@ -97,22 +96,19 @@ template<typename GeomTraits>
 using Col_4 = Eigen::Matrix<typename GeomTraits::FT, 4, 1>;
 
 template<typename VertexPointMap, typename TriangleMesh, typename GeomTraits>
-typename GeomTraits::Vector_3 construct_unit_normal_from_face(
+typename GeomTraits::Vector_3
-    const VertexPointMap& point_map,
+construct_unit_normal_from_face(typename boost::graph_traits<TriangleMesh>::face_descriptor f,
                                 const TriangleMesh& tmesh,
-    typename boost::graph_traits<TriangleMesh>::face_descriptor f, 
+                                const VertexPointMap point_map,
                                 const GeomTraits& gt)
 {
-  // initialize all necessary kernel functions
-  auto unit_normal = gt.construct_unit_normal_3_object();
-
-  // reference and descriptor types
   typedef typename boost::property_traits<VertexPointMap>::reference           point_reference;
   typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor      halfedge_descriptor;
 
+  auto unit_normal = gt.construct_unit_normal_3_object();
+
   const halfedge_descriptor h = halfedge(f, tmesh);
 
-  // get the three points of the face and calculate their unit normal
   const point_reference p = get(point_map, source(h, tmesh));
   const point_reference q = get(point_map, target(h, tmesh));
   const point_reference r = get(point_map, target(next(h, tmesh), tmesh));
@@ -120,17 +116,17 @@ typename GeomTraits::Vector_3 construct_unit_normal_from_face(
   return unit_normal(p, q, r);
 }
 
-template<typename VertexPointMap, typename TriangleMesh, typename GeomTraits> 
+template<typename TriangleMesh, typename VertexPointMap, typename GeomTraits>
-typename GeomTraits::Vector_3 construct_edge_normal(
+typename GeomTraits::Vector_3
-    const VertexPointMap& point_map,
+construct_edge_normal(typename boost::graph_traits<TriangleMesh>::halfedge_descriptor he,
                       const TriangleMesh& tmesh,
-    typename boost::graph_traits<TriangleMesh>::halfedge_descriptor he, 
+                      const VertexPointMap point_map,
                       const GeomTraits& gt)
 {
   typedef typename GeomTraits::Vector_3                                        Vector_3;
   typedef typename boost::graph_traits<TriangleMesh>::vertex_descriptor        vertex_descriptor;
 
-  const Vector_3 face_normal = construct_unit_normal_from_face(point_map, tmesh, face(he, tmesh), gt);
+  const Vector_3 face_normal = construct_unit_normal_from_face(face(he, tmesh), tmesh, point_map, gt);
 
   const vertex_descriptor vs = source(he, tmesh);
   const vertex_descriptor vt = target(he, tmesh);
@@ -139,15 +135,14 @@ typename GeomTraits::Vector_3 construct_edge_normal(
   const Vector_3 discontinuity_normal = cross_product(edge_vector, face_normal);
 
   // normalize
-  const Vector_3 normal = discontinuity_normal 
+  const Vector_3 normal = discontinuity_normal / sqrt(discontinuity_normal.squared_length());
-    / sqrt(discontinuity_normal.squared_length());
 
   return normal;
 }
 
 template<typename GeomTraits>
-Mat_4<GeomTraits> construct_classic_plane_quadric_from_normal(
+Mat_4<GeomTraits>
-    const typename GeomTraits::Vector_3& normal,
+construct_classic_plane_quadric_from_normal(const typename GeomTraits::Vector_3& normal,
                                             const typename GeomTraits::Point_3& point,
                                             const GeomTraits& gt)
 {
@@ -160,21 +155,20 @@ Mat_4<GeomTraits> construct_classic_plane_quadric_from_normal(
   const FT d = - dot_product(normal, construct_vector(ORIGIN, point));
 
   // row vector given by d appended to the normal
-  const Eigen::Matrix<FT, 1, 4> row (normal.x(), normal.y(), normal.z(), d);
+  const Eigen::Matrix<FT, 1, 4> row { normal.x(), normal.y(), normal.z(), d };
 
   // outer product
   return row.transpose() * row;
 }
 
-template<typename VertexPointMap, typename TriangleMesh, typename GeomTraits>
+template<typename TriangleMesh, typename VertexPointMap, typename GeomTraits>
-Mat_4<GeomTraits> construct_classic_plane_quadric_from_face(
+Mat_4<GeomTraits>
-    const VertexPointMap& point_map,
+construct_classic_plane_quadric_from_face(typename boost::graph_traits<TriangleMesh>::face_descriptor f,
                                           const TriangleMesh& mesh,
-    typename boost::graph_traits<TriangleMesh>::face_descriptor f,
+                                          const VertexPointMap point_map,
                                           const GeomTraits& gt)
 {
-  const typename GeomTraits::Vector_3 normal 
+  auto normal = construct_unit_normal_from_face(f, mesh, point_map, gt);
-    = construct_unit_normal_from_face(point_map, mesh, f, gt);
 
   // get any point of the face
   const auto p = get(point_map, source(halfedge(f, mesh), mesh));
@@ -182,18 +176,17 @@ Mat_4<GeomTraits> construct_classic_plane_quadric_from_face(
   return construct_classic_plane_quadric_from_normal(normal, p, gt);
 }
 
-template<typename VertexPointMap, typename TriangleMesh, typename GeomTraits>
+template<typename TriangleMesh, typename VertexPointMap, typename GeomTraits>
-Mat_4<GeomTraits> construct_classic_plane_quadric_from_edge(
+Mat_4<GeomTraits>
-    const VertexPointMap& point_map,
+construct_classic_plane_quadric_from_edge(typename boost::graph_traits<TriangleMesh>::halfedge_descriptor he,
                                           const TriangleMesh& mesh,
-    typename boost::graph_traits<TriangleMesh>::halfedge_descriptor he,
+                                          const VertexPointMap point_map,
                                           const GeomTraits& gt)
 {
-
   typedef typename GeomTraits::Vector_3                                        Vector_3;
   typedef typename boost::graph_traits<TriangleMesh>::vertex_descriptor        vertex_descriptor;
 
-  Vector_3 normal = construct_edge_normal(point_map, mesh, he, gt);
+  const Vector_3 normal = construct_edge_normal(he, mesh, point_map, gt);
 
   // use this normal to construct the quadric analogously to constructing quadric
   // from the normal of the face
@@ -201,21 +194,21 @@ Mat_4<GeomTraits> construct_classic_plane_quadric_from_edge(
 }
 
 template <typename GeomTraits>
-Mat_4<GeomTraits> construct_prob_plane_quadric_from_normal(
+Mat_4<GeomTraits>
-    const typename GeomTraits::Vector_3& mean_normal,
+construct_prob_plane_quadric_from_normal(const typename GeomTraits::Vector_3& mean_normal,
                                          const typename GeomTraits::Point_3& point,
                                          const GeomTraits& gt,
                                          typename GeomTraits::FT face_nv,
                                          typename GeomTraits::FT face_mv)
 {
-  auto squared_length = gt.compute_squared_length_3_object();
-  auto dot_product = gt.compute_scalar_product_3_object();
-  auto construct_vec_3 = gt.construct_vector_3_object();
-
   typedef typename GeomTraits::FT                                              FT;
   typedef typename GeomTraits::Vector_3                                        Vector_3;
   typedef Eigen::Matrix<FT, 4, 4>                                              Mat_4;
 
+  auto squared_length = gt.compute_squared_length_3_object();
+  auto dot_product = gt.compute_scalar_product_3_object();
+  auto construct_vec_3 = gt.construct_vector_3_object();
+
   const Vector_3 mean_vec = construct_vec_3(ORIGIN, point);
   const FT dot_mnmv = dot_product(mean_normal, mean_vec);
 
@@ -251,18 +244,16 @@ Mat_4<GeomTraits> construct_prob_plane_quadric_from_normal(
 }
 
 template<typename VertexPointMap, typename TriangleMesh, typename GeomTraits>
-std::array<typename GeomTraits::Vector_3, 3> vectors_from_face(
+std::array<typename GeomTraits::Vector_3, 3>
-    const VertexPointMap& point_map,
+vectors_from_face(const VertexPointMap point_map,
                   const TriangleMesh& tmesh,
                   typename boost::graph_traits<TriangleMesh>::face_descriptor f,
                   const GeomTraits& gt)
 {
-  auto construct_vector = gt.construct_vector_3_object();
-  
   typedef typename boost::property_traits<VertexPointMap>::reference           Point_reference;
   typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor      halfedge_descriptor;
 
-  std::array<typename GeomTraits::Vector_3, 3> arr { };
+  auto construct_vector = gt.construct_vector_3_object();
 
   const halfedge_descriptor h = halfedge(f, tmesh);
 
@@ -271,26 +262,26 @@ std::array<typename GeomTraits::Vector_3, 3> vectors_from_face(
   const Point_reference q = get(point_map, target(h, tmesh));
   const Point_reference r = get(point_map, target(next(h, tmesh), tmesh));
 
-  arr[0] = construct_vector(ORIGIN, p);
+  std::array<typename GeomTraits::Vector_3, 3> arr { construct_vector(ORIGIN, p),
-  arr[1] = construct_vector(ORIGIN, q);
+                                                     construct_vector(ORIGIN, q),
-  arr[2] = construct_vector(ORIGIN, r);
+                                                     construct_vector(ORIGIN, r) };
 
   return arr;
 }
 
-template<typename VertexPointMap, typename TriangleMesh, typename GeomTraits> 
+template<typename TriangleMesh, typename VertexPointMap, typename GeomTraits>
-Mat_4<GeomTraits> construct_classic_triangle_quadric_from_face(
+Mat_4<GeomTraits>
-    const VertexPointMap& point_map,
+construct_classic_triangle_quadric_from_face(typename boost::graph_traits<TriangleMesh>::face_descriptor f,
                                              const TriangleMesh& tmesh,
-    typename boost::graph_traits<TriangleMesh>::face_descriptor f, 
+                                             const VertexPointMap point_map,
                                              const GeomTraits& gt)
 {
+  typedef typename GeomTraits::FT                                              FT;
+
   auto cross_product = gt.construct_cross_product_vector_3_object();
   auto sum_vectors = gt.construct_sum_of_vectors_3_object();
   auto dot_product = gt.compute_scalar_product_3_object();
 
-  typedef typename GeomTraits::FT FT;
-
   auto vectors = vectors_from_face(point_map, tmesh, f, gt);
 
   Vector_3 a = vectors[0];
@@ -305,7 +296,6 @@ Mat_4<GeomTraits> construct_classic_triangle_quadric_from_face(
   const FT scalar_triple_product = dot_product(ab, c);
 
   Eigen::Matrix<FT, 1, 4> row;
-
   row << sum_of_cross_product.x(),   sum_of_cross_product.y(),
          sum_of_cross_product.z(), - scalar_triple_product;
 
@@ -314,8 +304,8 @@ Mat_4<GeomTraits> construct_classic_triangle_quadric_from_face(
 }
 
 template<typename GeomTraits>
-Eigen::Matrix<typename GeomTraits::FT, 3, 3> skew_sym_mat_cross_product(
+Eigen::Matrix<typename GeomTraits::FT, 3, 3>
-    const typename GeomTraits::Vector_3& v) 
+skew_sym_mat_cross_product(const typename GeomTraits::Vector_3& v)
 {
   Eigen::Matrix<typename GeomTraits::FT, 3, 3> mat;
 
@@ -326,14 +316,20 @@ Eigen::Matrix<typename GeomTraits::FT, 3, 3> skew_sym_mat_cross_product(
   return mat;
 }
 
-template<typename VertexPointMap, typename TriangleMesh, typename GeomTraits>
+template<typename TriangleMesh, typename VertexPointMap, typename GeomTraits>
-Mat_4<GeomTraits> construct_prob_triangle_quadric_from_face(
+Mat_4<GeomTraits>
-      const VertexPointMap& point_map,
+construct_prob_triangle_quadric_from_face(typename boost::graph_traits<TriangleMesh>::face_descriptor f,
-      const TriangleMesh& tmesh, 
-      typename boost::graph_traits<TriangleMesh>::face_descriptor f, 
                                           typename GeomTraits::FT var,
+                                          const TriangleMesh& tmesh,
+                                          const VertexPointMap point_map,
                                           const GeomTraits& gt)
 {
+  typedef typename GeomTraits::FT                                              FT;
+  typedef typename GeomTraits::Vector_3                                        Vector_3;
+
+  typedef Eigen::Matrix<FT, 3, 3>                                              Mat_3;
+  typedef Mat_4<GeomTraits>                                                    Mat_4;
+
   auto construct_vector = gt.construct_vector_3_object();
   auto cross_product = gt.construct_cross_product_vector_3_object();
   auto sum_vectors = gt.construct_sum_of_vectors_3_object();
@@ -343,14 +339,9 @@ Mat_4<GeomTraits> construct_prob_triangle_quadric_from_face(
   // the vertices of the given face
   auto vectors = vectors_from_face(point_map, tmesh, f, gt);
 
-  typedef typename GeomTraits::FT FT;
+  const Vector_3& a = vectors[0];
-  typedef typename GeomTraits::Vector_3 Vector_3;
+  const Vector_3& b = vectors[1];
-  typedef Eigen::Matrix<FT, 3, 3> Mat_3;
+  const Vector_3& c = vectors[2];
-  typedef Mat_4<GeomTraits> Mat_4;
-  
-  Vector_3 a = vectors[0];
-  Vector_3 b = vectors[1];
-  Vector_3 c = vectors[2];
 
   // calculate certain vectors used later
   const Vector_3 ab = cross_product(a, b);
@@ -374,33 +365,30 @@ Mat_4<GeomTraits> construct_prob_triangle_quadric_from_face(
   Mat_3 A = sum_cp_col * sum_cp_col.transpose();
   A += var * (cp_ab * cp_ab.transpose() + cp_bc * cp_bc.transpose() + cp_ca * cp_ca.transpose());
 
-  // add the 3 simple cross inference matrix - components (we only have one 
+  // add the 3 simple cross inference matrix - components (we only have one variance here)
-  // variance here)
+  A += 6 * square(var) * Mat_3::Identity();
-  A += 6 * var * var * Mat_3::Identity();
 
   // we need the determinant of matrix with columns a, b, c - we use the scalar triple product
   const FT det = dot_product(ab, c);
 
   // compute the b vector, this follows the formula directly - but we can factor
   // out the diagonal covariance matrices
-  const Eigen::Matrix<FT, 3, 1> res_b = det * sum_cp_col
+  const Eigen::Matrix<FT, 3, 1> res_b =
-    - var * (
+      det * sum_cp_col - var * (vector_to_col_vec<GeomTraits>(cross_product(a_minus_b, ab))
-        vector_to_col_vec<GeomTraits>(cross_product(a_minus_b, ab))
                                 + vector_to_col_vec<GeomTraits>(cross_product(b_minus_c, bc))
                                 + vector_to_col_vec<GeomTraits>(cross_product(c_minus_a, ca)))
-    + 2 * vector_to_col_vec<GeomTraits>(sum_vectors(sum_vectors(a, b), c)) * var * var;
+      + 2 * square(var) * vector_to_col_vec<GeomTraits>(sum_vectors(sum_vectors(a, b), c));
 
-  const FT res_c = det * det 
+  const FT ab2 = dot_product(ab, ab);
-    + var * (
+  const FT bc2 = dot_product(bc, bc);
-        dot_product(ab, ab) 
+  const FT ca2 = dot_product(ca, ca);
-        + dot_product(bc, bc) 
+  const FT a2 = dot_product(a, a);
-        + dot_product(ca, ca)) 
+  const FT b2 = dot_product(b, b);
-    + var * var * ( 
+  const FT c2 = dot_product(c, c);
-        2 * (
+
-          dot_product(a, a)
+  const FT res_c = square(det)
-          + dot_product(b, b)
+                   + var * (ab2 + bc2 + ca2)
-          + dot_product(c, c))
+                   + square(var) * (2 * (a2 + b2 + c2) + 6 * var);
-        + 6 * var);
 
   Mat_4 ret = Mat_4::Zero();
   ret.block(0, 0, 3, 3) = A;
@@ -424,8 +412,7 @@ Col_4<GeomTraits> construct_optimal_point_invertible(const Mat_4<GeomTraits>& qu
 }
 
 template <typename GeomTraits>
-Col_4<GeomTraits> construct_optimal_point_singular(
+Col_4<GeomTraits> construct_optimal_point_singular(const Mat_4<GeomTraits>& quadric,
-    const Mat_4<GeomTraits>& quadric,
                                                    const Col_4<GeomTraits>& p0,
                                                    const Col_4<GeomTraits>& p1)
 {
@@ -486,8 +473,9 @@ Col_4<GeomTraits> construct_optimal_point_singular(
 }
 
 template<typename TriangleMesh, typename GeomTraits>
-std::pair<typename GeomTraits::FT, typename GeomTraits::FT> estimate_variances(
+std::pair<typename GeomTraits::FT, typename GeomTraits::FT>
-    const TriangleMesh& mesh, const GeomTraits& gt, 
+estimate_variances(const TriangleMesh& mesh,
+                   const GeomTraits& gt,
                    typename GeomTraits::FT variance,
                    typename GeomTraits::FT p_factor)
 {

Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/GarlandHeckbert_policy_base.h

@@ -17,7 +17,7 @@
 
 #include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_profile.h>
 #include <CGAL/Surface_mesh_simplification/internal/Common.h>
-#include <iostream>
+
 #include <CGAL/tags.h>
 
 #include <Eigen/Dense>
@@ -29,22 +29,22 @@ namespace Surface_mesh_simplification {
 namespace internal {
 
 template <typename Mat_4>
-Mat_4 combine_matrices(const Mat_4& a, const Mat_4& b) {
+Mat_4 combine_matrices(const Mat_4& a, const Mat_4& b)
+{
   return a + b;
 }
 
-template<typename VCM,
+template <typename VertexCostMap,
           typename GeomTraits,
           typename QuadricImpl>
 class GarlandHeckbert_cost_base
 {
 public:
-
   // Tells the main function of 'Edge_collapse' that these
   // policies must call "initialize" and "update" functions.
   typedef CGAL::Tag_true                                                       Update_tag;
 
-  typedef VCM Vertex_cost_map;
+  typedef VertexCostMap                                                        Vertex_cost_map;
   typedef typename GeomTraits::FT                                              FT;
 
   typedef Eigen::Matrix<FT, 4, 4, Eigen::DontAlign>                            Mat_4;
@@ -53,27 +53,110 @@ class GarlandHeckbert_cost_base
   typedef typename GeomTraits::Point_3                                         Point_3;
   typedef typename GeomTraits::Vector_3                                        Vector_3;
 
-  GarlandHeckbert_cost_base() : m_cost_matrices() { }
+private:
+  Vertex_cost_map m_cost_matrices;
+
+  FT discontinuity_multiplier;
+
+public:
+  GarlandHeckbert_cost_base()
+    : m_cost_matrices()
+  { }
 
   GarlandHeckbert_cost_base(FT dm)
-    : m_cost_matrices(), discontinuity_multiplier(dm) { }
+    : m_cost_matrices(), discontinuity_multiplier(dm)
+  { }
 
   GarlandHeckbert_cost_base(Vertex_cost_map vcm, FT dm)
-    : m_cost_matrices(vcm), discontinuity_multiplier(dm) { }
+    : m_cost_matrices(vcm), discontinuity_multiplier(dm)
+  { }
+
+protected:
+  void init_vcm(const Vertex_cost_map vcm)
+  {
+    m_cost_matrices = vcm;
+  }
+
+  Col_4 point_to_homogenous_column(const Point_3& point) const
+  {
+    return Col_4(point.x(), point.y(), point.z(), FT(1));
+  }
+
+  template <typename TM, typename VPM>
+  Mat_4 construct_quadric(typename boost::graph_traits<TM>::face_descriptor f,
+                          const TM& tmesh,
+                          const VPM point_map,
+                          const GeomTraits& gt) const
+  {
+    return static_cast<const QuadricImpl*>(this)->construct_quadric_from_face(f, tmesh, point_map, gt);
+  }
+
+  template <typename TM, typename VPM>
+  Mat_4 construct_quadric(typename boost::graph_traits<TM>::halfedge_descriptor he,
+                          const TM& tmesh,
+                          const VPM point_map,
+                          const GeomTraits& gt) const
+  {
+    return static_cast<const QuadricImpl*>(this)->construct_quadric_from_edge(he, tmesh, point_map, gt);
+  }
+
+  template <typename TM, typename VPM>
+  Vector_3 construct_edge_normal(typename boost::graph_traits<TM>::halfedge_descriptor he,
+                                 const TM& tmesh,
+                                 const VPM point_map,
+                                 const GeomTraits& gt) const
+  {
+    typedef typename boost::graph_traits<TM>::vertex_descriptor vertex_descriptor;
+
+    const Vector_3 face_normal = construct_unit_normal_from_face(face(he, tmesh), tmesh, point_map, gt);
+
+    const vertex_descriptor vs = source(he, tmesh);
+    const vertex_descriptor vt = target(he, tmesh);
+
+    const Vector_3 edge_vector { get(point_map, vs), get(point_map, vt) };
+    const Vector_3 discontinuity_normal = cross_product(edge_vector, face_normal);
+
+    // normalize
+    const Vector_3 normal = discontinuity_normal / sqrt(discontinuity_normal.squared_length());
+
+    return normal;
+  }
+
+  template <typename VPM, typename TM>
+  Vector_3 construct_unit_normal_from_face(typename boost::graph_traits<TM>::face_descriptor f,
+                                           const TM& tmesh,
+                                           const VPM point_map,
+                                           const GeomTraits& gt) const
+  {
+    typedef typename boost::graph_traits<TM>::halfedge_descriptor              halfedge_descriptor;
+    typedef typename boost::property_traits<VPM>::reference                    Point_reference;
+
+    auto unit_normal = gt.construct_unit_normal_3_object();
+
+    const halfedge_descriptor h = halfedge(f, tmesh);
+
+    // get the three points of the face and calculate their unit normal
+    const Point_reference p = get(point_map, source(h, tmesh));
+    const Point_reference q = get(point_map, target(h, tmesh));
+    const Point_reference r = get(point_map, target(next(h, tmesh), tmesh));
+
+    return unit_normal(p, q, r);
+  }
 
   template <typename TM>
-  static bool is_discontinuity_edge(const typename boost::graph_traits<TM>::halfedge_descriptor& h,
+  static bool is_discontinuity_edge(const typename boost::graph_traits<TM>::halfedge_descriptor h,
                                     const TM& tmesh)
   {
     return is_border_edge(h, tmesh);
   }
 
+public:
   // initialize all quadrics
-  template<typename TM,
+  template <typename TM, typename VPM>
-           typename VPM>
+  void initialize(const TM& tmesh,
-  void initialize(const TM& tmesh, const VPM& vpm, const GeomTraits& gt) const
+                  const VPM vpm,
+                  const GeomTraits& gt) const
   {
-    // the graph library types used here
     typedef boost::graph_traits<TM>                                            GraphTraits;
     typedef typename GraphTraits::vertex_descriptor                            vertex_descriptor;
     typedef typename GraphTraits::halfedge_descriptor                          halfedge_descriptor;
@@ -83,21 +166,17 @@ class GarlandHeckbert_cost_base
     Mat_4 zero_mat = Mat_4::Zero();
 
     for(vertex_descriptor v : vertices(tmesh))
-    {
       put(m_cost_matrices, v, zero_mat);
-    }
 
     for(face_descriptor f : faces(tmesh))
     {
       if(f == GraphTraits::null_face())
-      {
         continue;
-      }
 
       const halfedge_descriptor h = halfedge(f, tmesh);
 
       // construtct the (4 x 4) matrix representing the plane quadric
-      const Mat_4 quadric = construct_quadric<VPM, TM>(vpm, tmesh, f, gt);
+      const Mat_4 quadric = construct_quadric(f, tmesh, vpm, gt);
 
       for(halfedge_descriptor shd : halfedges_around_face(h, tmesh))
       {
@@ -107,12 +186,10 @@ class GarlandHeckbert_cost_base
         put(m_cost_matrices, vs, combine_matrices(get(m_cost_matrices, vs), quadric));
 
         if(!is_discontinuity_edge(shd, tmesh))
-        {
           continue;
-        }
 
-        const Mat_4 discontinuous_quadric
+        const Mat_4 discontinuous_quadric =
-          = discontinuity_multiplier * construct_quadric<VPM, TM>(vpm, tmesh, shd, gt);
+          discontinuity_multiplier * construct_quadric(shd, tmesh, vpm, gt);
 
         put(m_cost_matrices, vs, combine_matrices(get(m_cost_matrices, vs), discontinuous_quadric));
         put(m_cost_matrices, vt, combine_matrices(get(m_cost_matrices, vt), discontinuous_quadric));
@@ -121,12 +198,14 @@ class GarlandHeckbert_cost_base
   }
 
   template <typename Profile>
-  boost::optional<typename Profile::FT> operator()(const Profile& profile,
+  boost::optional<typename Profile::FT>
+  operator()(const Profile& profile,
              const boost::optional<typename Profile::Point>& placement) const
   {
     typedef boost::optional<typename Profile::FT>                              Optional_FT;
 
-    if(!placement) {
+    if(!placement)
+    {
       // return empty value
       return boost::optional<typename Profile::FT>();
     }
@@ -150,101 +229,16 @@ class GarlandHeckbert_cost_base
         combine_matrices(get(m_cost_matrices, profile.v0()),
                          get(m_cost_matrices, profile.v1())));
   }
-
-  protected:
-  void init_vcm(const Vertex_cost_map& vcm) {
-    m_cost_matrices = vcm;
-  }
-  
-  template<typename VPM, typename TM>
-  Vector_3 construct_edge_normal(
-    const VPM& point_map,
-    const TM& tmesh, 
-    typename boost::graph_traits<TM>::halfedge_descriptor he, 
-    const GeomTraits& gt) const
-  {
-    typedef typename boost::graph_traits<TM>::vertex_descriptor vertex_descriptor;
-
-    const Vector_3 face_normal = this->construct_unit_normal_from_face(
-        point_map, tmesh, face(he, tmesh), gt);
-
-    const vertex_descriptor vs = source(he, tmesh);
-    const vertex_descriptor vt = target(he, tmesh);
-
-    const Vector_3 edge_vector = Vector_3(get(point_map, vs), get(point_map, vt));
-    const Vector_3 discontinuity_normal = cross_product(edge_vector, face_normal);
-
-    // normalize
-    const Vector_3 normal = discontinuity_normal 
-      / sqrt(discontinuity_normal.squared_length());
-
-    return normal;
-  }
-  
-  template<typename VPM, typename TM>
-    Vector_3 construct_unit_normal_from_face(
-        const VPM& point_map,
-        const TM& tmesh, 
-        typename boost::graph_traits<TM>::face_descriptor f, 
-        const GeomTraits& gt) const
-    {      
-      // initialize all necessary kernel functions
-      auto unit_normal = gt.construct_unit_normal_3_object();
-
-      // reference and descriptor types
-      typedef typename boost::property_traits<VPM>::reference Point_reference;
-      typedef typename boost::graph_traits<TM>::halfedge_descriptor halfedge_descriptor;
-
-      const halfedge_descriptor h = halfedge(f, tmesh);
-
-      // get the three points of the face and calculate their unit normal
-      const Point_reference p = get(point_map, source(h, tmesh));
-      const Point_reference q = get(point_map, target(h, tmesh));
-      const Point_reference r = get(point_map, target(next(h, tmesh), tmesh));
-
-      return unit_normal(p, q, r);
-    }
-
-  private:
-  Vertex_cost_map m_cost_matrices;
-
-  FT discontinuity_multiplier;
-
-  Col_4 point_to_homogenous_column(const Point_3& point) const 
-  {
-    return Col_4(point.x(), point.y(), point.z(), FT(1));
-  }
-
-  template<typename VPM, typename TM>
-  Mat_4 construct_quadric(
-      const VPM& point_map, 
-      const TM& tmesh,
-      typename boost::graph_traits<TM>::face_descriptor f, const GeomTraits& gt) const 
-  {
-    return static_cast<const QuadricImpl*>(this)
-      ->construct_quadric_from_face(point_map, tmesh, f, gt);
-  }
-
-  template<typename VPM, typename TM>
-  Mat_4 construct_quadric(
-      const VPM& point_map, 
-      const TM& tmesh,
-      typename boost::graph_traits<TM>::halfedge_descriptor he,
-      const GeomTraits& gt) const 
-  {
-    return static_cast<const QuadricImpl*>(this)
-      ->construct_quadric_from_edge(point_map, tmesh, he, gt);
-  }
 };
 
-template<typename VCM,
+template <typename VertexCostMap,
           typename GeomTraits,
           typename QuadricImpl>
 class GarlandHeckbert_placement_base
 {
 public:
   // define type required by the Get_cost concept
-    typedef VCM Vertex_cost_map;
+  typedef VertexCostMap                                                        Vertex_cost_map;
 
   // matrix and column vector types
   typedef typename GeomTraits::FT                                              FT;
@@ -252,32 +246,14 @@ template<typename VCM,
   typedef Eigen::Matrix<FT, 4, 4, Eigen::DontAlign>                            Mat_4;
   typedef Eigen::Matrix<FT, 4, 1>                                              Col_4;
 
+private:
+  Vertex_cost_map m_cost_matrices;
+
+public:
   GarlandHeckbert_placement_base() { }
-    GarlandHeckbert_placement_base(Vertex_cost_map cost_matrices) : m_cost_matrices(cost_matrices) { }
+  GarlandHeckbert_placement_base(Vertex_cost_map cost_matrices)
-
+    : m_cost_matrices(cost_matrices)
-    template<typename Profile>
+  { }
-    boost::optional<typename Profile::Point> operator()(const Profile& profile) const
-    {
-      CGAL_precondition(!get(m_cost_matrices, profile.v0()).isZero(0));
-      CGAL_precondition(!get(m_cost_matrices, profile.v1()).isZero(0));
-
-      // the combined matrix has already been computed in the evaluation of the cost...
-      const Mat_4 combinedMatrix = combine_matrices(
-          get(m_cost_matrices, profile.v0()),
-          get(m_cost_matrices, profile.v1()));
-
-      const Col_4 p0 = point_to_homogenous_column(profile.p0());
-      const Col_4 p1 = point_to_homogenous_column(profile.p1());
-
-      const Col_4 opt = construct_optimum(combinedMatrix, p0, p1);
-
-      boost::optional<typename Profile::Point> pt = typename Profile::Point(
-          opt(0) / opt(3),
-          opt(1) / opt(3),
-          opt(2) / opt(3));
-
-      return pt;
-    }
 
 protected:
   void init_vcm(const Vertex_cost_map& vcm)
@@ -285,10 +261,6 @@ template<typename VCM,
     m_cost_matrices = vcm;
   }
 
-
-  private:
-  Vertex_cost_map m_cost_matrices;
-
   // use CRTP to call the quadric implementation
   Col_4 construct_optimum(const Mat_4& mat, const Col_4& p0, const Col_4& p1) const
   {
@@ -299,6 +271,29 @@ template<typename VCM,
   {
     return Col_4(point.x(), point.y(), point.z(), FT(1));
   }
+
+public:
+  template <typename Profile>
+  boost::optional<typename Profile::Point> operator()(const Profile& profile) const
+  {
+    CGAL_precondition(!get(m_cost_matrices, profile.v0()).isZero(0));
+    CGAL_precondition(!get(m_cost_matrices, profile.v1()).isZero(0));
+
+    // the combined matrix has already been computed in the evaluation of the cost...
+    const Mat_4 combinedMatrix = combine_matrices(get(m_cost_matrices, profile.v0()),
+                                                  get(m_cost_matrices, profile.v1()));
+
+    const Col_4 p0 = point_to_homogenous_column(profile.p0());
+    const Col_4 p1 = point_to_homogenous_column(profile.p1());
+
+    const Col_4 opt = construct_optimum(combinedMatrix, p0, p1);
+
+    boost::optional<typename Profile::Point> pt = typename Profile::Point(opt(0) / opt(3),
+                                                                          opt(1) / opt(3),
+                                                                          opt(2) / opt(3));
+
+    return pt;
+  }
 };
 
 } // namespace internal