diff --git a/Surface_mesh_simplification/examples/Surface_mesh_simplification/gh_statistics.cpp b/Surface_mesh_simplification/examples/Surface_mesh_simplification/gh_statistics.cpp index d54056f1ff9..e95c76c1f31 100644 --- a/Surface_mesh_simplification/examples/Surface_mesh_simplification/gh_statistics.cpp +++ b/Surface_mesh_simplification/examples/Surface_mesh_simplification/gh_statistics.cpp @@ -90,9 +90,9 @@ void time_all_vector(const std::vector& meshes, const fs::path& ou // copying of the meshes std::vector temp_meshes = meshes; - for (int i = 0; i < n_burns; ++i) + for(int i = 0; i < n_burns; ++i) { - for (Surface_mesh& mesh : temp_meshes) + for(Surface_mesh& mesh : temp_meshes) { edge_collapse(mesh); } @@ -100,12 +100,12 @@ void time_all_vector(const std::vector& meshes, const fs::path& ou temp_meshes = meshes; // measure each run - for (int i = 0; i < n_samples; ++i) + for(int i = 0; i < n_samples; ++i) { // measure time taken by the edge_collapse function over each mesh individually time::time_point start_time = time::steady_clock::now(); - for (Surface_mesh& mesh : temp_meshes) + for(Surface_mesh& mesh : temp_meshes) { edge_collapse(mesh); } @@ -150,7 +150,7 @@ std::vector> get_all_meshes(const fs::path& Surface_mesh mesh; // iterate through all files in the given directory - for (fs::directory_iterator dir_iter(dir); dir_iter != end_iter; ++dir_iter) + for(fs::directory_iterator dir_iter(dir); dir_iter != end_iter; ++dir_iter) { // look for .off files if (fs::is_regular_file(dir_iter->status()) @@ -188,7 +188,7 @@ void time_policy(const fs::path& dir, const fs::path& output_file) auto vec = get_all_meshes(dir); std::vector meshes { }; - for (const auto& p : vec) { + for(const auto& p : vec) { std::cout << p.second << '\n'; } @@ -233,13 +233,13 @@ void hausdorff_errors_mesh(const Surface_mesh& mesh, const fs::path& out, InputI { fs::ofstream out_stream {out}; - for (InputIt it = begin; it != end; ++it) + for(InputIt it = begin; it != end; ++it) { std::array errs = hausdorff_errors(mesh, *it); out_stream << "ratio: " << *it << '\n'; - for (double e : errs) + for(double e : errs) { out_stream << std::to_string(e) << '\n'; } @@ -269,13 +269,13 @@ void hausdorff_errors_dir(const fs::path& dir, const fs::path& out, double ratio fs::ofstream out_stream {out}; - for (auto& p : meshes) + for(auto& p : meshes) { // calculate the errors and output to file std::array errors = hausdorff_errors(p.first, ratio); out_stream << p.second << '\n'; - for (double err : errors) + for(double err : errors) { out_stream << std::to_string(err) << '\n'; } @@ -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'; } diff --git a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_policies.h b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_policies.h index 45224a60753..ec39b18e978 100644 --- a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_policies.h +++ b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_policies.h @@ -24,97 +24,93 @@ 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 -class GarlandHeckbert_policies : +class GarlandHeckbert_policies : public internal::GarlandHeckbert_placement_base< - typename boost::property_map< - TriangleMesh, - CGAL::dynamic_vertex_property_t> - >::type, - GeomTraits, - GarlandHeckbert_policies - >, + typename boost::property_map< + TriangleMesh, + CGAL::dynamic_vertex_property_t > + >::type, + GeomTraits, + GarlandHeckbert_policies + >, public internal::GarlandHeckbert_cost_base< - typename boost::property_map< - TriangleMesh, - CGAL::dynamic_vertex_property_t> - >::type, - GeomTraits, - GarlandHeckbert_policies - > + typename boost::property_map< + TriangleMesh, + CGAL::dynamic_vertex_property_t > + >::type, + GeomTraits, + GarlandHeckbert_policies > { +public: + typedef typename GeomTraits::FT FT; - public: - typedef typename GeomTraits::FT FT; + // This is ugly, we later only use the Mat_4 from the Cost_base, but we want to define + // the matrix here already so it's nicer to define Cost_base in the first place + typedef typename Eigen::Matrix GH_matrix; + typedef CGAL::dynamic_vertex_property_t Cost_property; - // This is ugly, we later only use the Mat_4 from the - // Cost_base, but we want to define the matrix here already so it's nicer to define - // Cost_base in the first place - typedef typename Eigen::Matrix GH_matrix; - typedef CGAL::dynamic_vertex_property_t Cost_property; - - typedef typename boost::property_map::type Vertex_cost_map; + typedef typename boost::property_map::type Vertex_cost_map; typedef internal::GarlandHeckbert_placement_base< Vertex_cost_map, GeomTraits, GarlandHeckbert_policies - > Placement_base; + > Placement_base; typedef internal::GarlandHeckbert_cost_base< Vertex_cost_map, GeomTraits, GarlandHeckbert_policies - > Cost_base; + > Cost_base; - // both types are the same, this is so we avoid casting back to the base class in - // get_cost() or get_placement() - typedef GarlandHeckbert_policies Get_cost; - typedef GarlandHeckbert_policies Get_placement; + // both types are the same, this is so we avoid casting back to the base class in 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() - // 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 + // 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 - Mat_4 construct_quadric_from_face( - const VertexPointMap& point_map, - const TriangleMesh& tmesh, - typename boost::graph_traits::face_descriptor f, - const GeomTraits& gt) const + Mat_4 construct_quadric_from_face(typename boost::graph_traits::face_descriptor f, + const TriangleMesh& tmesh, + const VertexPointMap point_map, + const GeomTraits& gt) const { - return internal::construct_classic_plane_quadric_from_face( - point_map, tmesh, f, gt); + return internal::construct_classic_plane_quadric_from_face(f, tmesh, point_map, gt); } - + template - Mat_4 construct_quadric_from_edge( - const VertexPointMap& point_map, - const TriangleMesh& tmesh, - typename boost::graph_traits::halfedge_descriptor he, - const GeomTraits& gt) const + Mat_4 construct_quadric_from_edge(typename boost::graph_traits::halfedge_descriptor he, + const TriangleMesh& tmesh, + const VertexPointMap point_map, + const GeomTraits& gt) const { - return internal::construct_classic_plane_quadric_from_edge( - point_map, tmesh, he, gt); + return internal::construct_classic_plane_quadric_from_edge(he, tmesh, point_map, 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(quadric, p0, p1); } - const Get_cost& get_cost() const { return *this; } - const Get_placement& get_placement() const { return *this; } - - // 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 + 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(); @@ -151,7 +136,7 @@ class GarlandHeckbert_policies : } }; -} //namespace Surface_mesh_simplification -} //namespace CGAL +} // namespace Surface_mesh_simplification +} // namespace CGAL -#endif //CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_GARLANDHECKBERT_POLICIES_H +#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_GARLANDHECKBERT_POLICIES_H diff --git a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_probabilistic_policies.h b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_probabilistic_policies.h index ea15d9cac64..f91ecc40998 100644 --- a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_probabilistic_policies.h +++ b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_probabilistic_policies.h @@ -14,74 +14,63 @@ #define CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_GARLANDHECKBERT_PROBABILISTIC_POLICIES_H #include + #include #include #include + #include - + 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 probabilistic plane quadrics, optionally takes a face variance map -// giving a per-face variance -template::face_descriptor, - std::pair - >> -class GarlandHeckbert_probabilistic_policies : - public internal::GarlandHeckbert_placement_base< - typename boost::property_map< - TriangleMesh, - CGAL::dynamic_vertex_property_t> - >::type, - GeomTraits, - GarlandHeckbert_probabilistic_policies - >, - public internal::GarlandHeckbert_cost_base< - typename boost::property_map< - TriangleMesh, - CGAL::dynamic_vertex_property_t> - >::type, - GeomTraits, - GarlandHeckbert_probabilistic_policies - > +// Implements probabilistic plane quadrics, +// optionally takes a face variance map giving a per-face variance +template::face_descriptor, + std::pair > > +class GarlandHeckbert_probabilistic_policies + : public internal::GarlandHeckbert_placement_base< + typename boost::property_map< + TriangleMesh, + CGAL::dynamic_vertex_property_t > + >::type, + GeomTraits, + GarlandHeckbert_probabilistic_policies >, + public internal::GarlandHeckbert_cost_base< + typename boost::property_map< + TriangleMesh, + CGAL::dynamic_vertex_property_t> + >::type, + GeomTraits, + GarlandHeckbert_probabilistic_policies > { - - typedef typename GeomTraits::FT FT; - - typedef typename boost::graph_traits::face_descriptor face_descriptor; - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - - typedef typename boost::property_traits::value_type Face_variance; + typedef typename GeomTraits::FT FT; - public: + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + typedef typename boost::graph_traits::face_descriptor face_descriptor; - typedef typename Eigen::Matrix GH_matrix; - typedef CGAL::dynamic_vertex_property_t Cost_property; + typedef typename boost::property_traits::value_type Face_variance; - typedef typename boost::property_map::type Vertex_cost_map; +public: + typedef typename Eigen::Matrix GH_matrix; + typedef CGAL::dynamic_vertex_property_t Cost_property; + typedef typename boost::property_map::type Vertex_cost_map; typedef internal::GarlandHeckbert_placement_base< Vertex_cost_map, GeomTraits, GarlandHeckbert_probabilistic_policies - > Placement_base; + > Placement_base; typedef internal::GarlandHeckbert_cost_base< Vertex_cost_map, GeomTraits, GarlandHeckbert_probabilistic_policies - > Cost_base; + > Cost_base; // both types are the same, this is so we avoid casting back to the base class in // get_cost() or get_placement() - 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(); + typedef GarlandHeckbert_probabilistic_policies Get_cost; + typedef GarlandHeckbert_probabilistic_policies Get_placement; // 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) @@ -90,31 +79,45 @@ 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) + : GarlandHeckbert_probabilistic_policies(tmesh, 100) { } - + GarlandHeckbert_probabilistic_policies(TriangleMesh& tmesh, FT dm) : Cost_base(dm) - { + { // initialize the private variable vcm so it's lifetime is bound to that of the policy's vcm_ = get(Cost_property(), tmesh); // initialize both vcms Cost_base::init_vcm(vcm_); Placement_base::init_vcm(vcm_); - + // try to initialize the face variance map using the estimated variance // parameters are constants defined for this class - face_variance_map = FaceVarianceMap { - internal::estimate_variances(tmesh, GeomTraits(), good_default_variance_unit, - position_variance_factor) }; + face_variance_map = + FaceVarianceMap { internal::estimate_variances(tmesh, GeomTraits(), + good_default_variance_unit, + position_variance_factor) }; } - + GarlandHeckbert_probabilistic_policies(TriangleMesh& tmesh, FT dm, const FaceVarianceMap& fvm) : Cost_base(dm), face_variance_map(fvm) - { + { // initialize the private variable vcm so it's lifetime is bound to that of the policy's vcm_ = get(Cost_property(), tmesh); @@ -122,70 +125,56 @@ class GarlandHeckbert_probabilistic_policies : 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; } + + // 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 - Mat_4 construct_quadric_from_face( - const VPM& point_map, - const TriangleMesh& tmesh, - face_descriptor f, - const GeomTraits& gt) const + Mat_4 construct_quadric_from_face(face_descriptor f, + const TriangleMesh& tmesh, + const VPM point_map, + const GeomTraits& gt) const { - const Vector_3 normal = internal::construct_unit_normal_from_face< - VPM, TriangleMesh, GeomTraits>(point_map, tmesh, f, gt); - + const Vector_3 normal = internal::construct_unit_normal_from_face(f, tmesh, point_map, 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 - Mat_4 construct_quadric_from_edge( - const VPM& point_map, - const TriangleMesh& tmesh, - halfedge_descriptor he, - const GeomTraits& gt) const - { - const Vector_3 normal = internal::construct_edge_normal(point_map, tmesh, he, gt); + template + Mat_4 construct_quadric_from_edge(halfedge_descriptor he, + const TriangleMesh& tmesh, + const VPM point_map, + const GeomTraits& gt) const + { + 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, - const Col_4& p1) const + + Col_4 construct_optimal_point(const Mat_4& quadric, + const Col_4& p0, + const Col_4& p1) const { return internal::construct_optimal_point_invertible(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 } // namespace CGAL -#endif //CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_GARLANDHECKBERT_PROBABILISTIC_POLICIES_H +#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_GARLANDHECKBERT_PROBABILISTIC_POLICIES_H diff --git a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_probabilistic_tri_policies.h b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_probabilistic_tri_policies.h index 427cd085717..c867388632d 100644 --- a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_probabilistic_tri_policies.h +++ b/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 + #include #include + #include #include @@ -28,57 +30,50 @@ 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::face_descriptor, typename GeomTraits::FT - >> -class GarlandHeckbert_probabilistic_tri_policies : - public internal::GarlandHeckbert_placement_base< - typename boost::property_map< - TriangleMesh, - CGAL::dynamic_vertex_property_t> - >::type, - GeomTraits, - GarlandHeckbert_probabilistic_tri_policies - >, +template::face_descriptor, + typename GeomTraits::FT> > +class GarlandHeckbert_probabilistic_tri_policies + : public internal::GarlandHeckbert_placement_base< + typename boost::property_map< + TriangleMesh, + CGAL::dynamic_vertex_property_t> + >::type, + GeomTraits, + GarlandHeckbert_probabilistic_tri_policies + >, public internal::GarlandHeckbert_cost_base< - typename boost::property_map< - TriangleMesh, - CGAL::dynamic_vertex_property_t> - >::type, - GeomTraits, - GarlandHeckbert_probabilistic_tri_policies - > + typename boost::property_map< + TriangleMesh, + CGAL::dynamic_vertex_property_t> + >::type, + GeomTraits, + GarlandHeckbert_probabilistic_tri_policies > { + typedef typename boost::property_traits::value_type Face_variance; - typedef typename boost::property_traits::value_type Face_variance; +public: + typedef typename GeomTraits::FT FT; - public: - typedef typename GeomTraits::FT FT; - - typedef typename Eigen::Matrix Mat_3; - typedef typename Eigen::Matrix GH_matrix; - typedef CGAL::dynamic_vertex_property_t Cost_property; - - typedef typename boost::property_map::type Vertex_cost_map; + typedef typename Eigen::Matrix Mat_3; + typedef typename Eigen::Matrix GH_matrix; + typedef CGAL::dynamic_vertex_property_t Cost_property; + typedef typename boost::property_map::type Vertex_cost_map; typedef internal::GarlandHeckbert_placement_base< Vertex_cost_map, GeomTraits, GarlandHeckbert_probabilistic_tri_policies - > Placement_base; + > Placement_base; typedef internal::GarlandHeckbert_cost_base< Vertex_cost_map, GeomTraits, GarlandHeckbert_probabilistic_tri_policies - > Cost_base; + > Cost_base; // both types are the same, this is so we avoid casting back to the base class in // get_cost() or get_placement() - 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(); + typedef GarlandHeckbert_probabilistic_tri_policies Get_cost; + typedef GarlandHeckbert_probabilistic_tri_policies Get_placement; // 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) @@ -86,15 +81,27 @@ class GarlandHeckbert_probabilistic_tri_policies : using typename Cost_base::Col_4; 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) + : GarlandHeckbert_probabilistic_tri_policies(tmesh, 100) { } - + GarlandHeckbert_probabilistic_tri_policies(TriangleMesh& tmesh, FT dm) : Cost_base(dm) - { + { // initialize the private variable vcm so it's lifetime is bound to that of the policy's vcm_ = get(Cost_property(), tmesh); @@ -107,19 +114,19 @@ 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 }; } GarlandHeckbert_probabilistic_tri_policies(TriangleMesh& tmesh, - FT dm, - const FaceVarianceMap* fvm) + FT dm, + 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); @@ -127,61 +134,51 @@ class GarlandHeckbert_probabilistic_tri_policies : 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; } + + // 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 - Mat_4 construct_quadric_from_face( - const VPM& point_map, - const TM& tmesh, - typename boost::graph_traits::face_descriptor f, - const GeomTraits& gt) const + Mat_4 construct_quadric_from_face(typename boost::graph_traits::face_descriptor f, + const TM& tmesh, + const VPM point_map, + const GeomTraits& gt) const { FT variance = get(face_variance_map, f); - return internal::construct_prob_triangle_quadric_from_face( - point_map, tmesh, f, variance, gt); + return internal::construct_prob_triangle_quadric_from_face(f, variance, tmesh, point_map, 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 - Mat_4 construct_quadric_from_edge( - const VPM& point_map, - const TM& tmesh, - typename boost::graph_traits::halfedge_descriptor he, - const GeomTraits& gt) const + Mat_4 construct_quadric_from_edge(typename boost::graph_traits::halfedge_descriptor he, + const TM& tmesh, + 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)); - - return internal::construct_prob_plane_quadric_from_normal(normal, p, gt, variance, - position_variance_factor * variance); + + return internal::construct_prob_plane_quadric_from_normal(normal, p, gt, variance, + position_variance_factor * variance); } - + Col_4 construct_optimal_point(const Mat_4 quadric, const Col_4& p0, const Col_4& p1) const { return internal::construct_optimal_point_invertible(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 } // namespace CGAL -#endif //CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_GARLANDHECKBERT_PROBABILISTIC_TRI_POLICIES_H +#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_GARLANDHECKBERT_PROBABILISTIC_TRI_POLICIES_H diff --git a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_triangle_policies.h b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_triangle_policies.h index 787977bdf8c..02fc626ca09 100644 --- a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_triangle_policies.h +++ b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_triangle_policies.h @@ -18,55 +18,60 @@ #include #include #include + #include + #include namespace CGAL { namespace Surface_mesh_simplification { -// use triangle quadrics for the faces, classical plane quadrics for the edges +// use triangle quadrics for the faces, classical plane quadrics for the edges // and optimize with a check for singular matrices // // implements classic triangle policies template -class GarlandHeckbert_triangle_policies : - public internal::GarlandHeckbert_cost_base< - typename boost::property_map< - TriangleMesh, - CGAL::dynamic_vertex_property_t> - >::type, - GeomTraits, - GarlandHeckbert_triangle_policies - >, +class GarlandHeckbert_triangle_policies + : public internal::GarlandHeckbert_cost_base< + typename boost::property_map< + TriangleMesh, + CGAL::dynamic_vertex_property_t> + >::type, + GeomTraits, + GarlandHeckbert_triangle_policies + >, public internal::GarlandHeckbert_placement_base< - typename boost::property_map< - TriangleMesh, - CGAL::dynamic_vertex_property_t> - >::type, - GeomTraits, - GarlandHeckbert_triangle_policies - > + typename boost::property_map< + TriangleMesh, + CGAL::dynamic_vertex_property_t> + >::type, + GeomTraits, + GarlandHeckbert_triangle_policies > { - public: - typedef typename GeomTraits::FT FT; - - typedef Eigen::Matrix Mat_4; - typedef Eigen::Matrix Col_4; - - typedef CGAL::dynamic_vertex_property_t Cost_property; - typedef typename boost::property_map::type Vertex_cost_map; - +public: + typedef typename GeomTraits::FT FT; + + typedef Eigen::Matrix Mat_4; + typedef Eigen::Matrix Col_4; + + typedef CGAL::dynamic_vertex_property_t Cost_property; + typedef typename boost::property_map::type Vertex_cost_map; + typedef internal::GarlandHeckbert_placement_base< Vertex_cost_map, GeomTraits, GarlandHeckbert_triangle_policies - > Placement_base; + > Placement_base; typedef internal::GarlandHeckbert_cost_base< Vertex_cost_map, GeomTraits, GarlandHeckbert_triangle_policies - > Cost_base; - - using Cost_base::operator(); - using Placement_base::operator(); - + > Cost_base; + + typedef GarlandHeckbert_triangle_policies Get_cost; + 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,45 +80,39 @@ 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 - Mat_4 construct_quadric_from_face( - const VertexPointMap& point_map, - const TriangleMesh& tmesh, - typename boost::graph_traits::face_descriptor f, - const GeomTraits& gt) const + Mat_4 construct_quadric_from_face(typename boost::graph_traits::face_descriptor f, + const TriangleMesh& tmesh, + const VertexPointMap point_map, + const GeomTraits& gt) const { - return internal::construct_classic_triangle_quadric_from_face( - point_map, tmesh, f, gt); + return internal::construct_classic_triangle_quadric_from_face(f, tmesh, point_map, gt); } - + template - Mat_4 construct_quadric_from_edge( - const VertexPointMap& point_map, - const TriangleMesh& tmesh, - typename boost::graph_traits::halfedge_descriptor he, - const GeomTraits& gt) const + Mat_4 construct_quadric_from_edge(typename boost::graph_traits::halfedge_descriptor he, + const TriangleMesh& tmesh, + const VertexPointMap point_map, + const GeomTraits& gt) const { - return internal::construct_classic_plane_quadric_from_edge( - point_map, tmesh, he, gt); + return internal::construct_classic_plane_quadric_from_edge(he, tmesh, point_map, 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(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 -} //namespace CGAL +} // namespace Surface_mesh_simplification +} // namespace CGAL -#endif //CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_GARLANDHECKBERT_POLICIES_H +#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_GARLANDHECKBERT_POLICIES_H diff --git a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/GarlandHeckbert_functions.h b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/GarlandHeckbert_functions.h index 7487bc35e07..05878cfc2f1 100644 --- a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/GarlandHeckbert_functions.h +++ b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/GarlandHeckbert_functions.h @@ -17,6 +17,7 @@ #include #include + #include namespace CGAL { @@ -28,7 +29,7 @@ template bool invert_matrix_4(const Matrix& m, Matrix& im) { double det; - + double A2323 = m(2, 2) * m(3, 3) - m(2, 3) * m(3, 2); double A1323 = m(2, 1) * m(3, 3) - m(2, 3) * m(3, 1); double A1223 = m(2, 1) * m(3, 2) - m(2, 2) * m(3, 1); @@ -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 -Eigen::Matrix vector_to_col_vec( - const typename GeomTraits::Vector_3& v) +Eigen::Matrix +vector_to_col_vec(const typename GeomTraits::Vector_3& v) { Eigen::Matrix col {v.x(), v.y(), v.z()}; return col; @@ -97,22 +96,19 @@ template using Col_4 = Eigen::Matrix; template -typename GeomTraits::Vector_3 construct_unit_normal_from_face( - const VertexPointMap& point_map, - const TriangleMesh& tmesh, - typename boost::graph_traits::face_descriptor f, - const GeomTraits& gt) -{ - // initialize all necessary kernel functions - auto unit_normal = gt.construct_unit_normal_3_object(); +typename GeomTraits::Vector_3 +construct_unit_normal_from_face(typename boost::graph_traits::face_descriptor f, + const TriangleMesh& tmesh, + const VertexPointMap point_map, + const GeomTraits& gt) +{ + typedef typename boost::property_traits::reference point_reference; + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - // reference and descriptor types - typedef typename boost::property_traits::reference point_reference; - typedef typename boost::graph_traits::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 GeomTraits::Vector_3 construct_edge_normal( - const VertexPointMap& point_map, - const TriangleMesh& tmesh, - typename boost::graph_traits::halfedge_descriptor he, - const GeomTraits& gt) +template +typename GeomTraits::Vector_3 +construct_edge_normal(typename boost::graph_traits::halfedge_descriptor he, + const TriangleMesh& tmesh, + const VertexPointMap point_map, + const GeomTraits& gt) { - typedef typename GeomTraits::Vector_3 Vector_3; - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename GeomTraits::Vector_3 Vector_3; + typedef typename boost::graph_traits::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,17 +135,16 @@ 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 - / sqrt(discontinuity_normal.squared_length()); + const Vector_3 normal = discontinuity_normal / sqrt(discontinuity_normal.squared_length()); return normal; } template -Mat_4 construct_classic_plane_quadric_from_normal( - const typename GeomTraits::Vector_3& normal, - const typename GeomTraits::Point_3& point, - const GeomTraits& gt) +Mat_4 +construct_classic_plane_quadric_from_normal(const typename GeomTraits::Vector_3& normal, + const typename GeomTraits::Point_3& point, + const GeomTraits& gt) { typedef typename GeomTraits::FT FT; @@ -157,70 +152,68 @@ Mat_4 construct_classic_plane_quadric_from_normal( auto construct_vector = gt.construct_vector_3_object(); // negative dot product between the normal and the position vector - const FT d = -dot_product(normal, construct_vector(ORIGIN, point)); + const FT d = - dot_product(normal, construct_vector(ORIGIN, point)); // row vector given by d appended to the normal - const Eigen::Matrix row (normal.x(), normal.y(), normal.z(), d); + const Eigen::Matrix row { normal.x(), normal.y(), normal.z(), d }; // outer product return row.transpose() * row; } -template -Mat_4 construct_classic_plane_quadric_from_face( - const VertexPointMap& point_map, - const TriangleMesh& mesh, - typename boost::graph_traits::face_descriptor f, - const GeomTraits& gt) +template +Mat_4 +construct_classic_plane_quadric_from_face(typename boost::graph_traits::face_descriptor f, + const TriangleMesh& mesh, + const VertexPointMap point_map, + const GeomTraits& gt) { - const typename GeomTraits::Vector_3 normal - = construct_unit_normal_from_face(point_map, mesh, f, gt); + auto normal = construct_unit_normal_from_face(f, mesh, point_map, gt); - // get any point of the face + // get any point of the face const auto p = get(point_map, source(halfedge(f, mesh), mesh)); return construct_classic_plane_quadric_from_normal(normal, p, gt); } -template -Mat_4 construct_classic_plane_quadric_from_edge( - const VertexPointMap& point_map, - const TriangleMesh& mesh, - typename boost::graph_traits::halfedge_descriptor he, - const GeomTraits& gt) +template +Mat_4 +construct_classic_plane_quadric_from_edge(typename boost::graph_traits::halfedge_descriptor he, + const TriangleMesh& mesh, + const VertexPointMap point_map, + const GeomTraits& gt) { + typedef typename GeomTraits::Vector_3 Vector_3; + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - typedef typename GeomTraits::Vector_3 Vector_3; - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + const Vector_3 normal = construct_edge_normal(he, mesh, point_map, gt); - Vector_3 normal = construct_edge_normal(point_map, mesh, he, gt); - // use this normal to construct the quadric analogously to constructing quadric // from the normal of the face return construct_classic_plane_quadric_from_normal(normal, get(point_map, source(he, mesh)), gt); } template -Mat_4 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) +Mat_4 +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) { + typedef typename GeomTraits::FT FT; + typedef typename GeomTraits::Vector_3 Vector_3; + typedef Eigen::Matrix 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(); - typedef typename GeomTraits::FT FT; - typedef typename GeomTraits::Vector_3 Vector_3; - typedef Eigen::Matrix Mat_4; - const Vector_3 mean_vec = construct_vec_3(ORIGIN, point); const FT dot_mnmv = dot_product(mean_normal, mean_vec); // Eigen column vector of length 3 - const Eigen::Matrix mean_n_col {mean_normal.x(), mean_normal.y(), mean_normal.z()}; + const Eigen::Matrix mean_n_col { mean_normal.x(), mean_normal.y(), mean_normal.z() }; // start by setting values along the diagonal Mat_4 mat = face_nv * Mat_4::Identity(); @@ -235,7 +228,7 @@ Mat_4 construct_prob_plane_quadric_from_normal( // the b column and row appear with a negative sign const auto b1 = -(dot_mnmv * mean_normal + face_nv * mean_vec); - const Eigen::Matrix b {b1.x(), b1.y(), b1.z()}; + const Eigen::Matrix b { b1.x(), b1.y(), b1.z() }; mat.col(3).head(3) = b; mat.row(3).head(3) = b.transpose(); @@ -243,27 +236,25 @@ Mat_4 construct_prob_plane_quadric_from_normal( // set the value in the bottom right corner, we get this by considering // that we only have single variances given instead of covariance matrices mat(3, 3) = CGAL::square(dot_mnmv) - + face_nv * squared_length(mean_vec) - + face_mv * squared_length(mean_normal) - + 3 * face_nv * face_mv; + + face_nv * squared_length(mean_vec) + + face_mv * squared_length(mean_normal) + + 3 * face_nv * face_mv; return mat; } template -std::array vectors_from_face( - const VertexPointMap& point_map, - const TriangleMesh& tmesh, - typename boost::graph_traits::face_descriptor f, - const GeomTraits& gt) +std::array +vectors_from_face(const VertexPointMap point_map, + const TriangleMesh& tmesh, + typename boost::graph_traits::face_descriptor f, + const GeomTraits& gt) { + typedef typename boost::property_traits::reference Point_reference; + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + auto construct_vector = gt.construct_vector_3_object(); - - typedef typename boost::property_traits::reference Point_reference; - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - - std::array arr { }; - + const halfedge_descriptor h = halfedge(f, tmesh); // get all points and turn them into location vectors so we can use cross product on them @@ -271,28 +262,28 @@ std::array 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); - arr[1] = construct_vector(ORIGIN, q); - arr[2] = construct_vector(ORIGIN, r); + std::array arr { construct_vector(ORIGIN, p), + construct_vector(ORIGIN, q), + construct_vector(ORIGIN, r) }; - return arr; + return arr; } -template -Mat_4 construct_classic_triangle_quadric_from_face( - const VertexPointMap& point_map, - const TriangleMesh& tmesh, - typename boost::graph_traits::face_descriptor f, - const GeomTraits& gt) +template +Mat_4 +construct_classic_triangle_quadric_from_face(typename boost::graph_traits::face_descriptor f, + const TriangleMesh& tmesh, + 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]; Vector_3 b = vectors[1]; Vector_3 c = vectors[2]; @@ -305,52 +296,52 @@ Mat_4 construct_classic_triangle_quadric_from_face( const FT scalar_triple_product = dot_product(ab, c); Eigen::Matrix row; - - row << sum_of_cross_product.x(), sum_of_cross_product.y(), - sum_of_cross_product.z(), -scalar_triple_product; + row << sum_of_cross_product.x(), sum_of_cross_product.y(), + sum_of_cross_product.z(), - scalar_triple_product; // calculate the outer product of row^t*row return row.transpose() * row; } template -Eigen::Matrix skew_sym_mat_cross_product( - const typename GeomTraits::Vector_3& v) +Eigen::Matrix +skew_sym_mat_cross_product(const typename GeomTraits::Vector_3& v) { Eigen::Matrix mat; - - mat << 0, -v.z(), v.y(), - v.z(), 0, -v.x(), - -v.y(), v.x(), 0; + + mat << 0, - v.z(), v.y(), + v.z(), 0, - v.x(), + - v.y(), v.x(), 0; return mat; } -template -Mat_4 construct_prob_triangle_quadric_from_face( - const VertexPointMap& point_map, - const TriangleMesh& tmesh, - typename boost::graph_traits::face_descriptor f, - typename GeomTraits::FT var, - const GeomTraits& gt) +template +Mat_4 +construct_prob_triangle_quadric_from_face(typename boost::graph_traits::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 Mat_3; + typedef Mat_4 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(); auto dot_product = gt.compute_scalar_product_3_object(); - + // array containing the position vectors corresponding to // the vertices of the given face auto vectors = vectors_from_face(point_map, tmesh, f, gt); - - typedef typename GeomTraits::FT FT; - typedef typename GeomTraits::Vector_3 Vector_3; - typedef Eigen::Matrix Mat_3; - typedef Mat_4 Mat_4; - - Vector_3 a = vectors[0]; - Vector_3 b = vectors[1]; - Vector_3 c = vectors[2]; + + const Vector_3& a = vectors[0]; + const Vector_3& b = vectors[1]; + const Vector_3& c = vectors[2]; // calculate certain vectors used later const Vector_3 ab = cross_product(a, b); @@ -368,44 +359,41 @@ Mat_4 construct_prob_triangle_quadric_from_face( const Vector_3 sum_of_cross_product = sum_vectors(sum_vectors(ab, bc), ca); - const Eigen::Matrix - sum_cp_col{ sum_of_cross_product.x(), sum_of_cross_product.y(), sum_of_cross_product.z() }; + const Eigen::Matrix + sum_cp_col { sum_of_cross_product.x(), sum_of_cross_product.y(), sum_of_cross_product.z() }; 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 - // variance here) - A += 6 * var * var * Mat_3::Identity(); + // add the 3 simple cross inference matrix - components (we only have one variance here) + A += 6 * square(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 res_b = det * sum_cp_col - - var * ( - vector_to_col_vec(cross_product(a_minus_b, ab)) - + vector_to_col_vec(cross_product(b_minus_c, bc)) - + vector_to_col_vec(cross_product(c_minus_a, ca))) - + 2 * vector_to_col_vec(sum_vectors(sum_vectors(a, b), c)) * var * var; + // out the diagonal covariance matrices + const Eigen::Matrix res_b = + det * sum_cp_col - var * (vector_to_col_vec(cross_product(a_minus_b, ab)) + + vector_to_col_vec(cross_product(b_minus_c, bc)) + + vector_to_col_vec(cross_product(c_minus_a, ca))) + + 2 * square(var) * vector_to_col_vec(sum_vectors(sum_vectors(a, b), c)); - const FT res_c = det * det - + var * ( - dot_product(ab, ab) - + dot_product(bc, bc) - + dot_product(ca, ca)) - + var * var * ( - 2 * ( - dot_product(a, a) - + dot_product(b, b) - + dot_product(c, c)) - + 6 * var); + const FT ab2 = dot_product(ab, ab); + const FT bc2 = dot_product(bc, bc); + const FT ca2 = dot_product(ca, ca); + const FT a2 = dot_product(a, a); + const FT b2 = dot_product(b, b); + const FT c2 = dot_product(c, c); + + const FT res_c = square(det) + + var * (ab2 + bc2 + ca2) + + square(var) * (2 * (a2 + b2 + c2) + 6 * var); Mat_4 ret = Mat_4::Zero(); ret.block(0, 0, 3, 3) = A; - ret.block(3, 0, 1, 3) = -res_b.transpose(); - ret.block(0, 3, 3, 1) = -res_b; + ret.block(3, 0, 1, 3) = - res_b.transpose(); + ret.block(0, 3, 3, 1) = - res_b; ret(3, 3) = res_c; return ret; @@ -424,13 +412,12 @@ Col_4 construct_optimal_point_invertible(const Mat_4& qu } template -Col_4 construct_optimal_point_singular( - const Mat_4& quadric, - const Col_4& p0, - const Col_4& p1) +Col_4 construct_optimal_point_singular(const Mat_4& quadric, + const Col_4& p0, + const Col_4& p1) { - typedef typename GeomTraits::FT FT; - + typedef typename GeomTraits::FT FT; + // in this case, the matrix mat may no be invertible, // so we save the result to check Mat_4 mat; @@ -438,12 +425,12 @@ Col_4 construct_optimal_point_singular( Mat_4 inverse; bool invertible = invert_matrix_4(mat, inverse); - - if (invertible) + + if(invertible) { return inverse.col(3); } - else + else { Col_4 opt_pt; @@ -451,11 +438,11 @@ Col_4 construct_optimal_point_singular( const FT a = (p1mp0.transpose() * quadric * p1mp0)(0, 0); const FT b = 2 * (p0.transpose() * quadric * p1mp0)(0, 0); - if (a == 0) + if(a == 0) { - if (b < 0) + if(b < 0) opt_pt = p1; - else if (b == 0) + else if(b == 0) opt_pt = 0.5 * (p0 + p1); else opt_pt = p0; @@ -486,27 +473,28 @@ Col_4 construct_optimal_point_singular( } template -std::pair estimate_variances( - const TriangleMesh& mesh, const GeomTraits& gt, - typename GeomTraits::FT variance, - typename GeomTraits::FT p_factor) +std::pair +estimate_variances(const TriangleMesh& mesh, + const GeomTraits& gt, + typename GeomTraits::FT variance, + typename GeomTraits::FT p_factor) { - typedef typename TriangleMesh::Vertex_index vertex_descriptor; - typedef typename TriangleMesh::Edge_index edge_descriptor; - typedef typename GeomTraits::FT FT; - typedef typename GeomTraits::Point_3 Point_3; + typedef typename TriangleMesh::Vertex_index vertex_descriptor; + typedef typename TriangleMesh::Edge_index edge_descriptor; + typedef typename GeomTraits::FT FT; + typedef typename GeomTraits::Point_3 Point_3; FT average_edge_length = 0; auto construct_vector = gt.construct_vector_3_object(); - for (edge_descriptor e : edges(mesh)) + for(edge_descriptor e : edges(mesh)) { vertex_descriptor v1 = mesh.vertex(e, 0); vertex_descriptor v2 = mesh.vertex(e, 1); - const Point_3& p1 = mesh.point(v1); - const Point_3& p2 = mesh.point(v2); + const Point_3& p1 = mesh.point(v1); + const Point_3& p2 = mesh.point(v2); const Vector_3 vec = construct_vector(p1, p2); average_edge_length += sqrt(vec.squared_length()); diff --git a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/GarlandHeckbert_policy_base.h b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/GarlandHeckbert_policy_base.h index e59092ce276..dfe6156d28c 100644 --- a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/GarlandHeckbert_policy_base.h +++ b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/GarlandHeckbert_policy_base.h @@ -17,7 +17,7 @@ #include #include -#include + #include #include @@ -28,76 +28,155 @@ namespace CGAL { namespace Surface_mesh_simplification { namespace internal { -template -Mat_4 combine_matrices(const Mat_4& a, const Mat_4& b) { +template +Mat_4 combine_matrices(const Mat_4& a, const Mat_4& b) +{ return a + b; } -template +template class GarlandHeckbert_cost_base { - public: - +public: // Tells the main function of 'Edge_collapse' that these // policies must call "initialize" and "update" functions. - typedef CGAL::Tag_true Update_tag; + typedef CGAL::Tag_true Update_tag; - typedef VCM Vertex_cost_map; - typedef typename GeomTraits::FT FT; + typedef VertexCostMap Vertex_cost_map; + typedef typename GeomTraits::FT FT; - typedef Eigen::Matrix Mat_4; - typedef Eigen::Matrix Col_4; + typedef Eigen::Matrix Mat_4; + typedef Eigen::Matrix Col_4; - typedef typename GeomTraits::Point_3 Point_3; - typedef typename GeomTraits::Vector_3 Vector_3; + 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) + { } - template - static bool is_discontinuity_edge(const typename boost::graph_traits::halfedge_descriptor& h, +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 + Mat_4 construct_quadric(typename boost::graph_traits::face_descriptor f, + const TM& tmesh, + const VPM point_map, + const GeomTraits& gt) const + { + return static_cast(this)->construct_quadric_from_face(f, tmesh, point_map, gt); + } + + template + Mat_4 construct_quadric(typename boost::graph_traits::halfedge_descriptor he, + const TM& tmesh, + const VPM point_map, + const GeomTraits& gt) const + { + return static_cast(this)->construct_quadric_from_edge(he, tmesh, point_map, gt); + } + + template + Vector_3 construct_edge_normal(typename boost::graph_traits::halfedge_descriptor he, + const TM& tmesh, + const VPM point_map, + const GeomTraits& gt) const + { + typedef typename boost::graph_traits::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 + Vector_3 construct_unit_normal_from_face(typename boost::graph_traits::face_descriptor f, + const TM& tmesh, + const VPM point_map, + const GeomTraits& gt) const + { + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + typedef typename boost::property_traits::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 + static bool is_discontinuity_edge(const typename boost::graph_traits::halfedge_descriptor h, const TM& tmesh) { return is_border_edge(h, tmesh); } +public: // initialize all quadrics - template - void initialize(const TM& tmesh, const VPM& vpm, const GeomTraits& gt) const + template + void initialize(const TM& tmesh, + const VPM vpm, + const GeomTraits& gt) const { - // the graph library types used here - typedef boost::graph_traits GraphTraits; - typedef typename GraphTraits::vertex_descriptor vertex_descriptor; - typedef typename GraphTraits::halfedge_descriptor halfedge_descriptor; - typedef typename GraphTraits::face_descriptor face_descriptor; - typedef typename boost::property_traits::reference Point_reference; + typedef boost::graph_traits GraphTraits; + typedef typename GraphTraits::vertex_descriptor vertex_descriptor; + typedef typename GraphTraits::halfedge_descriptor halfedge_descriptor; + typedef typename GraphTraits::face_descriptor face_descriptor; + typedef typename boost::property_traits::reference Point_reference; 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()) - { + 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, 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 - = discontinuity_multiplier * construct_quadric(vpm, tmesh, shd, gt); + const Mat_4 discontinuous_quadric = + 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)); @@ -120,13 +197,15 @@ class GarlandHeckbert_cost_base } } - template - boost::optional operator()(const Profile& profile, - const boost::optional& placement) const + template + boost::optional + operator()(const Profile& profile, + const boost::optional& placement) const { - typedef boost::optional Optional_FT; + typedef boost::optional Optional_FT; - if(!placement) { + if(!placement) + { // return empty value return boost::optional(); } @@ -142,163 +221,79 @@ class GarlandHeckbert_cost_base return cost; } - template + template void update_after_collapse(const Profile& profile, - const VertexDescriptor new_v) const + const VertexDescriptor new_v) const { put(m_cost_matrices, new_v, 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 - Vector_3 construct_edge_normal( - const VPM& point_map, - const TM& tmesh, - typename boost::graph_traits::halfedge_descriptor he, - const GeomTraits& gt) const - { - typedef typename boost::graph_traits::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 - Vector_3 construct_unit_normal_from_face( - const VPM& point_map, - const TM& tmesh, - typename boost::graph_traits::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::reference Point_reference; - typedef typename boost::graph_traits::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 - Mat_4 construct_quadric( - const VPM& point_map, - const TM& tmesh, - typename boost::graph_traits::face_descriptor f, const GeomTraits& gt) const - { - return static_cast(this) - ->construct_quadric_from_face(point_map, tmesh, f, gt); - } - - template - Mat_4 construct_quadric( - const VPM& point_map, - const TM& tmesh, - typename boost::graph_traits::halfedge_descriptor he, - const GeomTraits& gt) const - { - return static_cast(this) - ->construct_quadric_from_edge(point_map, tmesh, he, gt); - } }; -template - class GarlandHeckbert_placement_base +template +class GarlandHeckbert_placement_base { - public: - // define type required by the Get_cost concept - typedef VCM Vertex_cost_map; +public: + // define type required by the Get_cost concept + typedef VertexCostMap Vertex_cost_map; - // matrix and column vector types - typedef typename GeomTraits::FT FT; - typedef typename GeomTraits::Point_3 Point_3; - typedef Eigen::Matrix Mat_4; - typedef Eigen::Matrix Col_4; + // matrix and column vector types + typedef typename GeomTraits::FT FT; + typedef typename GeomTraits::Point_3 Point_3; + typedef Eigen::Matrix Mat_4; + typedef Eigen::Matrix Col_4; - GarlandHeckbert_placement_base() { } - GarlandHeckbert_placement_base(Vertex_cost_map cost_matrices) : m_cost_matrices(cost_matrices) { } +private: + Vertex_cost_map m_cost_matrices; - template - boost::optional 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)); +public: + GarlandHeckbert_placement_base() { } + GarlandHeckbert_placement_base(Vertex_cost_map cost_matrices) + : m_cost_matrices(cost_matrices) + { } - // 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 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) +protected: + void init_vcm(const Vertex_cost_map& 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 + Col_4 construct_optimum(const Mat_4& mat, const Col_4& p0, const Col_4& p1) const { return static_cast(this)->construct_optimal_point(mat, p0, p1); } - Col_4 point_to_homogenous_column(const Point_3& point) const + Col_4 point_to_homogenous_column(const Point_3& point) const { return Col_4(point.x(), point.y(), point.z(), FT(1)); } + +public: + template + boost::optional 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 pt = typename Profile::Point(opt(0) / opt(3), + opt(1) / opt(3), + opt(2) / opt(3)); + + return pt; + } }; } // namespace internal