From b342e119674d403571b940ba339997ac6ef76e7b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Wed, 22 Dec 2021 12:03:11 +0100 Subject: [PATCH] Refactor + improvements for Garland Heckbert methods --- .../Edge_collapse/GarlandHeckbert_policies.h | 153 ++-- .../GarlandHeckbert_probabilistic_policies.h | 212 +++--- ...rlandHeckbert_probabilistic_tri_policies.h | 220 +++--- .../GarlandHeckbert_triangle_policies.h | 134 ++-- .../internal/GarlandHeckbert_functions.h | 691 ++++++++++-------- .../internal/GarlandHeckbert_policy_base.h | 283 ++++--- 6 files changed, 848 insertions(+), 845 deletions(-) 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 ec39b18e978..4782502f4ca 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 @@ -16,95 +16,22 @@ #include #include -#include #include - -#include +#include namespace CGAL { namespace Surface_mesh_simplification { -template -class GarlandHeckbert_policies : - public internal::GarlandHeckbert_placement_base< - 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 > +template +class Plane_quadric_calculator { -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; - - typedef typename boost::property_map::type Vertex_cost_map; - - typedef internal::GarlandHeckbert_placement_base< - Vertex_cost_map, GeomTraits, GarlandHeckbert_policies - > Placement_base; - - typedef internal::GarlandHeckbert_cost_base< - Vertex_cost_map, GeomTraits, GarlandHeckbert_policies - > 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; - - // 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; - -private: - Vertex_cost_map vcm_; + typedef typename internal::GarlandHeckbert_matrix_types::Mat_4 Mat_4; + typedef typename internal::GarlandHeckbert_matrix_types::Col_4 Col_4; public: - GarlandHeckbert_policies(TriangleMesh& tmesh, - FT dm = FT(100)) // unused @tmp - { - vcm_ = get(Cost_property(), tmesh); + Plane_quadric_calculator() { } - // 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(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(f, tmesh, point_map, gt); - } - - template + template Mat_4 construct_quadric_from_edge(typename boost::graph_traits::halfedge_descriptor he, const TriangleMesh& tmesh, const VertexPointMap point_map, @@ -113,15 +40,22 @@ public: 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 + template + 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_optimal_point_singular(quadric, p0, p1); + return internal::construct_classic_plane_quadric_from_face(f, tmesh, point_map, gt); } - Mat_4 construct_quadric_from_normal(const Vector_3& normal, - const Point_3& point, + // @fixme unused? + Mat_4 construct_quadric_from_normal(const typename GeomTraits::Vector_3& normal, + const typename GeomTraits::Point_3& point, const GeomTraits& gt) const { + typedef typename GeomTraits::FT FT; + auto dot_product = gt.compute_scalar_product_3_object(); auto construct_vector = gt.construct_vector_3_object(); @@ -129,11 +63,60 @@ public: 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 Col_4 row { normal.x(), normal.y(), normal.z(), d }; // outer product return row.transpose() * row; } + + 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); + } +}; + +template +class GarlandHeckbert_policies + : public internal::GarlandHeckbert_placement_base< + Plane_quadric_calculator, TriangleMesh, GeomTraits>, + public internal::GarlandHeckbert_cost_base< + Plane_quadric_calculator, TriangleMesh, GeomTraits> +{ +public: + typedef Plane_quadric_calculator Quadric_calculator; + +private: + typedef internal::GarlandHeckbert_placement_base< + Quadric_calculator, TriangleMesh, GeomTraits> Placement_base; + typedef internal::GarlandHeckbert_cost_base< + Quadric_calculator, TriangleMesh, GeomTraits> Cost_base; + + // Diamond base + typedef internal::GarlandHeckbert_quadrics_storage< + Quadric_calculator, TriangleMesh, GeomTraits> Quadrics_storage; + + typedef GarlandHeckbert_policies Self; + +public: + typedef Self Get_cost; + typedef Self Get_placement; + + typedef typename GeomTraits::FT FT; + +public: + GarlandHeckbert_policies(TriangleMesh& tmesh, + const FT dm = FT(100)) + : Quadrics_storage(tmesh), Placement_base(), Cost_base(dm) + { } + +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(); }; } // namespace Surface_mesh_simplification 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 f91ecc40998..225c7889e74 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 @@ -19,161 +19,165 @@ #include #include -#include +#include +#include + +#include +#include namespace CGAL { namespace Surface_mesh_simplification { -// 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 > +template +class Probabilistic_plane_quadric_calculator { - typedef typename GeomTraits::FT FT; - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; typedef typename boost::graph_traits::face_descriptor face_descriptor; - typedef typename boost::property_traits::value_type Face_variance; + typedef typename GeomTraits::FT FT; -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 Constant_property_map > Default_FVM; + typedef typename Default::Get::type Face_variance_map; - typedef internal::GarlandHeckbert_placement_base< - Vertex_cost_map, GeomTraits, GarlandHeckbert_probabilistic_policies - > Placement_base; - - typedef internal::GarlandHeckbert_cost_base< - Vertex_cost_map, GeomTraits, GarlandHeckbert_probabilistic_policies - > 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; - - // 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; + typedef typename internal::GarlandHeckbert_matrix_types::Mat_4 Mat_4; + typedef typename internal::GarlandHeckbert_matrix_types::Col_4 Col_4; private: - Vertex_cost_map vcm_; + // @fixme check the magic values // magic number determined by some testing - static constexpr FT good_default_variance_unit = 0.05; + static constexpr FT 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; + Face_variance_map m_face_variance_map; public: - // default discontinuity multiplier is 100 - GarlandHeckbert_probabilistic_policies(TriangleMesh& tmesh) - : GarlandHeckbert_probabilistic_policies(tmesh, 100) + Probabilistic_plane_quadric_calculator() { CGAL_assertion(false); } // compilation only + + template + Probabilistic_plane_quadric_calculator(const FVM fvm) + : m_face_variance_map(fvm) { } - GarlandHeckbert_probabilistic_policies(TriangleMesh& tmesh, FT dm) - : Cost_base(dm) + Probabilistic_plane_quadric_calculator(TriangleMesh& tmesh, + typename boost::enable_if >::type* = nullptr) { - // 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) }; + m_face_variance_map = Default_FVM { internal::estimate_variances(tmesh, GeomTraits(), + 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); - - // initialize both vcms - 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(face_descriptor f, + template + Mat_4 construct_quadric_from_edge(const halfedge_descriptor he, const TriangleMesh& tmesh, - const VPM point_map, + const VertexPointMap vpm, const GeomTraits& gt) const { - 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)); + typedef typename GeomTraits::FT FT; + typedef typename GeomTraits::Point_3 Point_3; + typedef typename GeomTraits::Vector_3 Vector_3; - FT n_variance; - FT p_variance; - std::tie(n_variance, p_variance) = get(face_variance_map, f); + const Vector_3 normal = internal::construct_edge_normal(he, tmesh, vpm, gt); + const Point_3 p = get(vpm, source(he, tmesh)); + + FT n_variance, p_variance; + std::tie(n_variance, p_variance) = get(m_face_variance_map, face(he, tmesh)); return internal::construct_prob_plane_quadric_from_normal(normal, p, gt, n_variance, p_variance); } - template - Mat_4 construct_quadric_from_edge(halfedge_descriptor he, + template + Mat_4 construct_quadric_from_face(const face_descriptor f, const TriangleMesh& tmesh, - const VPM point_map, + const VertexPointMap vpm, 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)); + typedef typename GeomTraits::FT FT; + typedef typename GeomTraits::Point_3 Point_3; + typedef typename GeomTraits::Vector_3 Vector_3; - FT n_variance; - FT p_variance; - std::tie(n_variance, p_variance) = get(face_variance_map, face(he, tmesh)); + const Vector_3 normal = internal::construct_unit_normal_from_face(f, tmesh, vpm, gt); + const Point_3 p = get(vpm, source(halfedge(f, tmesh), tmesh)); + + FT n_variance, p_variance; + std::tie(n_variance, p_variance) = get(m_face_variance_map, f); 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 + const Col_4& /*p0*/, + const Col_4& /*p1*/) const { + // @fixme check this return internal::construct_optimal_point_invertible(quadric); } }; +// Implements probabilistic plane quadrics, +// optionally takes a face variance map giving a per-face variance +template +class GarlandHeckbert_probabilistic_policies + : public internal::GarlandHeckbert_placement_base< + Probabilistic_plane_quadric_calculator, + TriangleMesh, GeomTraits>, + public internal::GarlandHeckbert_cost_base< + Probabilistic_plane_quadric_calculator, + TriangleMesh, GeomTraits> +{ +public: + typedef Probabilistic_plane_quadric_calculator< + TriangleMesh, GeomTraits, FaceVarianceMap> Quadric_calculator; + +private: + typedef internal::GarlandHeckbert_placement_base< + Quadric_calculator, TriangleMesh, GeomTraits> Placement_base; + + typedef internal::GarlandHeckbert_cost_base< + Quadric_calculator, TriangleMesh, GeomTraits> Cost_base; + + // Diamond base + typedef internal::GarlandHeckbert_quadrics_storage< + Quadric_calculator, TriangleMesh, GeomTraits> Quadrics_storage; + + typedef GarlandHeckbert_probabilistic_policies Self; + +public: + typedef Self Get_cost; + typedef Self Get_placement; + + typedef typename GeomTraits::FT FT; + +public: + // Only available if the quadric calculator is using the default (constant) variance property map + GarlandHeckbert_probabilistic_policies(TriangleMesh& tmesh, + const FT dm = FT(100)) + : Quadrics_storage(tmesh, Quadric_calculator(tmesh)), Placement_base(), Cost_base(dm) + { } + + template + GarlandHeckbert_probabilistic_policies(TriangleMesh& tmesh, + const FT dm, + const FVM fvm) + : Quadrics_storage(tmesh, Quadric_calculator(fvm)), Placement_base(), Cost_base(dm) + { } + +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(); +}; + } // namespace Surface_mesh_simplification } // namespace CGAL 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 c867388632d..bcc950b9f19 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 @@ -17,167 +17,155 @@ #include #include -#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" -// -// 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 - >, - public internal::GarlandHeckbert_cost_base< - typename boost::property_map< - TriangleMesh, - CGAL::dynamic_vertex_property_t> - >::type, - GeomTraits, - GarlandHeckbert_probabilistic_tri_policies > +template +class Probabilistic_triangle_quadric_calculator { - typedef typename boost::property_traits::value_type Face_variance; - -public: typedef typename GeomTraits::FT FT; + typedef typename GeomTraits::Point_3 Point_3; + typedef typename GeomTraits::Vector_3 Vector_3; - 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 boost::graph_traits::face_descriptor face_descriptor; - typedef internal::GarlandHeckbert_placement_base< - Vertex_cost_map, GeomTraits, GarlandHeckbert_probabilistic_tri_policies - > Placement_base; + typedef Constant_property_map Default_FVM; + typedef typename Default::Get::type Face_variance_map; - typedef internal::GarlandHeckbert_cost_base< - Vertex_cost_map, GeomTraits, GarlandHeckbert_probabilistic_tri_policies - > 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; - - // 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; + typedef typename internal::GarlandHeckbert_matrix_types::Mat_4 Mat_4; + typedef typename internal::GarlandHeckbert_matrix_types::Col_4 Col_4; 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; + static constexpr FT 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; + Face_variance_map m_face_variance_map; + public: - // default discontinuity multiplier is 100 - GarlandHeckbert_probabilistic_tri_policies(TriangleMesh& tmesh) - : GarlandHeckbert_probabilistic_tri_policies(tmesh, 100) + Probabilistic_triangle_quadric_calculator() { CGAL_assertion(false); } // compilation only + + template + Probabilistic_triangle_quadric_calculator(const FVM fvm) + : m_face_variance_map(fvm) { } - GarlandHeckbert_probabilistic_tri_policies(TriangleMesh& tmesh, FT dm) - : Cost_base(dm) + Probabilistic_triangle_quadric_calculator(TriangleMesh& tmesh, + typename boost::enable_if >::type* = nullptr) { - // 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_); - - FT variance; - FT discard_position; - - std::tie(variance, discard_position) = internal::estimate_variances(tmesh, GeomTraits(), - good_default_variance_unit, position_variance_factor); + // try to initialize the face variance map using the estimated variance + // parameters are constants defined for this class + FT variance, discard_position; + std::tie(variance, discard_position) = + internal::estimate_variances(tmesh, GeomTraits(), default_variance_unit, position_variance_factor); // see probabilistic plane quadrics - face_variance_map = FaceVarianceMap { variance }; + m_face_variance_map = Default_FVM { variance }; } - GarlandHeckbert_probabilistic_tri_policies(TriangleMesh& tmesh, - 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. && 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); - - // initialize both vcms - 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(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(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(typename boost::graph_traits::halfedge_descriptor he, - const TM& tmesh, - const VPM point_map, + template + Mat_4 construct_quadric_from_edge(typename boost::graph_traits::halfedge_descriptor he, + const TriangleMesh& tmesh, + const VertexPointMap point_map, const GeomTraits& gt) const { // same as in probabilistic plane policy 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)); + const FT variance = get(m_face_variance_map, face(he, tmesh)); + // @fixme plane? 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 + template + Mat_4 construct_quadric_from_face(typename boost::graph_traits::face_descriptor f, + const TriangleMesh& tmesh, + const VertexPointMap point_map, + const GeomTraits& gt) const { + const FT variance = get(m_face_variance_map, f); + + return internal::construct_prob_triangle_quadric_from_face(f, variance, tmesh, point_map, gt); + } + + Col_4 construct_optimal_point(const Mat_4& quadric, const Col_4& /*p0*/, const Col_4& /*p1*/) const + { + // @fixme check this return internal::construct_optimal_point_invertible(quadric); } }; +// implements probabilistic triangle faces and optionally takes a face variance map +// analogously to probabilistic plane quadrics +template +class GarlandHeckbert_probabilistic_tri_policies + : public internal::GarlandHeckbert_placement_base< + Probabilistic_triangle_quadric_calculator, + TriangleMesh, GeomTraits>, + public internal::GarlandHeckbert_cost_base< + Probabilistic_triangle_quadric_calculator, + TriangleMesh, GeomTraits> +{ +public: + typedef Probabilistic_triangle_quadric_calculator< + TriangleMesh, GeomTraits, FaceVarianceMap> Quadric_calculator; + +private: + typedef internal::GarlandHeckbert_placement_base< + Quadric_calculator, TriangleMesh, GeomTraits> Placement_base; + + typedef internal::GarlandHeckbert_cost_base< + Quadric_calculator, TriangleMesh, GeomTraits> Cost_base; + + typedef internal::GarlandHeckbert_quadrics_storage< + Quadric_calculator, TriangleMesh, GeomTraits> Quadrics_storage; + + typedef GarlandHeckbert_probabilistic_tri_policies Self; + +public: + typedef Self Get_cost; + typedef Self Get_placement; + + typedef typename GeomTraits::FT FT; + +public: + // Only available if the quadric calculator is using the default (constant) variance property map + GarlandHeckbert_probabilistic_tri_policies(TriangleMesh& tmesh, + const FT dm = FT(100)) + : Quadrics_storage(tmesh, Quadric_calculator(tmesh)), Placement_base(), Cost_base(dm) + { } + + template + GarlandHeckbert_probabilistic_tri_policies(TriangleMesh& tmesh, + const FT dm, + const FVM fvm) + : Quadrics_storage(tmesh, Quadric_calculator(fvm)), Placement_base(), Cost_base(dm) + { } + +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(); +}; + } // namespace Surface_mesh_simplification } // namespace CGAL 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 02fc626ca09..b578c63e841 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 @@ -19,76 +19,29 @@ #include #include -#include - -#include - namespace CGAL { namespace Surface_mesh_simplification { -// 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 - >, - public internal::GarlandHeckbert_placement_base< - typename boost::property_map< - TriangleMesh, - CGAL::dynamic_vertex_property_t> - >::type, - GeomTraits, - GarlandHeckbert_triangle_policies > +template +class Triangle_quadric_calculator { -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; - - typedef internal::GarlandHeckbert_cost_base< - Vertex_cost_map, GeomTraits, GarlandHeckbert_triangle_policies - > Cost_base; - - typedef GarlandHeckbert_triangle_policies Get_cost; - typedef GarlandHeckbert_triangle_policies Get_placement; - -private: - Vertex_cost_map vcm_; + typedef typename internal::GarlandHeckbert_matrix_types::Mat_4 Mat_4; + typedef typename internal::GarlandHeckbert_matrix_types::Col_4 Col_4; public: - GarlandHeckbert_triangle_policies(TriangleMesh& tmesh, FT dm = FT(100)) + Triangle_quadric_calculator() { } + + template + Mat_4 construct_quadric_from_edge(typename boost::graph_traits::halfedge_descriptor he, + const TriangleMesh& tmesh, + const VertexPointMap point_map, + const GeomTraits& gt) const { - vcm_ = get(Cost_property(), tmesh); - - Cost_base::init_vcm(vcm_); - Placement_base::init_vcm(vcm_); + // @fixme "plane"? why not incident triangles? + return internal::construct_classic_plane_quadric_from_edge(he, tmesh, point_map, gt); } -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 + template Mat_4 construct_quadric_from_face(typename boost::graph_traits::face_descriptor f, const TriangleMesh& tmesh, const VertexPointMap point_map, @@ -97,21 +50,60 @@ public: return internal::construct_classic_triangle_quadric_from_face(f, tmesh, point_map, gt); } - template - 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(he, tmesh, point_map, gt); - } - - 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_singular(quadric, p0, p1); } }; +// 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_placement_base< + Triangle_quadric_calculator, TriangleMesh, GeomTraits>, + public internal::GarlandHeckbert_cost_base< + Triangle_quadric_calculator, TriangleMesh, GeomTraits> +{ +public: + typedef Triangle_quadric_calculator Quadric_calculator; + +private: + typedef internal::GarlandHeckbert_placement_base< + Quadric_calculator, TriangleMesh, GeomTraits> Placement_base; + typedef internal::GarlandHeckbert_cost_base< + Quadric_calculator, TriangleMesh, GeomTraits> Cost_base; + + // Diamond base + typedef internal::GarlandHeckbert_quadrics_storage< + Quadric_calculator, TriangleMesh, GeomTraits> Quadrics_storage; + + typedef GarlandHeckbert_triangle_policies Self; + +public: + typedef Self Get_cost; + typedef Self Get_placement; + + typedef typename GeomTraits::FT FT; + +public: + GarlandHeckbert_triangle_policies(TriangleMesh& tmesh, + const FT dm = FT(100)) + : Quadrics_storage(tmesh), Placement_base(), Cost_base(dm) + { } + +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(); +}; + } // namespace Surface_mesh_simplification } // namespace CGAL 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 05878cfc2f1..ab0049a4f4e 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 @@ -16,7 +16,7 @@ #include -#include +#include #include @@ -25,289 +25,166 @@ namespace Surface_mesh_simplification { namespace internal { // taken from https://stackoverflow.com/questions/1148309/inverting-a-4x4-matrix/ -template +template bool invert_matrix_4(const Matrix& m, Matrix& im) { - double det; + typedef typename Matrix::value_type FT; - 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); - double A0323 = m(2, 0) * m(3, 3) - m(2, 3) * m(3, 0); - double A0223 = m(2, 0) * m(3, 2) - m(2, 2) * m(3, 0); - double A0123 = m(2, 0) * m(3, 1) - m(2, 1) * m(3, 0); - double A2313 = m(1, 2) * m(3, 3) - m(1, 3) * m(3, 2); - double A1313 = m(1, 1) * m(3, 3) - m(1, 3) * m(3, 1); - double A1213 = m(1, 1) * m(3, 2) - m(1, 2) * m(3, 1); - double A2312 = m(1, 2) * m(2, 3) - m(1, 3) * m(2, 2); - double A1312 = m(1, 1) * m(2, 3) - m(1, 3) * m(2, 1); - double A1212 = m(1, 1) * m(2, 2) - m(1, 2) * m(2, 1); - double A0313 = m(1, 0) * m(3, 3) - m(1, 3) * m(3, 0); - double A0213 = m(1, 0) * m(3, 2) - m(1, 2) * m(3, 0); - double A0312 = m(1, 0) * m(2, 3) - m(1, 3) * m(2, 0); - double A0212 = m(1, 0) * m(2, 2) - m(1, 2) * m(2, 0); - double A0113 = m(1, 0) * m(3, 1) - m(1, 1) * m(3, 0); - double A0112 = m(1, 0) * m(2, 1) - m(1, 1) * m(2, 0); + FT det; + + FT A2323 = m(2, 2) * m(3, 3) - m(2, 3) * m(3, 2); + FT A1323 = m(2, 1) * m(3, 3) - m(2, 3) * m(3, 1); + FT A1223 = m(2, 1) * m(3, 2) - m(2, 2) * m(3, 1); + FT A0323 = m(2, 0) * m(3, 3) - m(2, 3) * m(3, 0); + FT A0223 = m(2, 0) * m(3, 2) - m(2, 2) * m(3, 0); + FT A0123 = m(2, 0) * m(3, 1) - m(2, 1) * m(3, 0); + // FT A2313 = m(1, 2) * m(3, 3) - m(1, 3) * m(3, 2); + // FT A1313 = m(1, 1) * m(3, 3) - m(1, 3) * m(3, 1); + // FT A1213 = m(1, 1) * m(3, 2) - m(1, 2) * m(3, 1); + FT A2312 = m(1, 2) * m(2, 3) - m(1, 3) * m(2, 2); + FT A1312 = m(1, 1) * m(2, 3) - m(1, 3) * m(2, 1); + FT A1212 = m(1, 1) * m(2, 2) - m(1, 2) * m(2, 1); + // FT A0313 = m(1, 0) * m(3, 3) - m(1, 3) * m(3, 0); + // FT A0213 = m(1, 0) * m(3, 2) - m(1, 2) * m(3, 0); + FT A0312 = m(1, 0) * m(2, 3) - m(1, 3) * m(2, 0); + FT A0212 = m(1, 0) * m(2, 2) - m(1, 2) * m(2, 0); + // FT A0113 = m(1, 0) * m(3, 1) - m(1, 1) * m(3, 0); + FT A0112 = m(1, 0) * m(2, 1) - m(1, 1) * m(2, 0); det = m(0, 0) * ( m(1, 1) * A2323 - m(1, 2) * A1323 + m(1, 3) * A1223 ) - m(0, 1) * ( m(1, 0) * A2323 - m(1, 2) * A0323 + m(1, 3) * A0223 ) + m(0, 2) * ( m(1, 0) * A1323 - m(1, 1) * A0323 + m(1, 3) * A0123 ) - m(0, 3) * ( m(1, 0) * A1223 - m(1, 1) * A0223 + m(1, 2) * A0123 ); - det = 1 / det; - if(det == 0.) + if(is_zero(det)) return false; + det = 1 / det; + // we never actually use values other than those in the third column, // so might as well not calculate them - //im(0, 0) = det * ( m(1, 1) * A2323 - m(1, 2) * A1323 + m(1, 3) * A1223 ); - //im(0, 1) = det * - ( m(0, 1) * A2323 - m(0, 2) * A1323 + m(0, 3) * A1223 ); - //im(0, 2) = det * ( m(0, 1) * A2313 - m(0, 2) * A1313 + m(0, 3) * A1213 ); + // im(0, 0) = det * ( m(1, 1) * A2323 - m(1, 2) * A1323 + m(1, 3) * A1223 ); + // im(0, 1) = det * - ( m(0, 1) * A2323 - m(0, 2) * A1323 + m(0, 3) * A1223 ); + // im(0, 2) = det * ( m(0, 1) * A2313 - m(0, 2) * A1313 + m(0, 3) * A1213 ); im(0, 3) = det * - ( m(0, 1) * A2312 - m(0, 2) * A1312 + m(0, 3) * A1212 ); - //im(1, 0) = det * - ( m(1, 0) * A2323 - m(1, 2) * A0323 + m(1, 3) * A0223 ); - //im(1, 1) = det * ( m(0, 0) * A2323 - m(0, 2) * A0323 + m(0, 3) * A0223 ); - //im(1, 2) = det * - ( m(0, 0) * A2313 - m(0, 2) * A0313 + m(0, 3) * A0213 ); + // im(1, 0) = det * - ( m(1, 0) * A2323 - m(1, 2) * A0323 + m(1, 3) * A0223 ); + // im(1, 1) = det * ( m(0, 0) * A2323 - m(0, 2) * A0323 + m(0, 3) * A0223 ); + // im(1, 2) = det * - ( m(0, 0) * A2313 - m(0, 2) * A0313 + m(0, 3) * A0213 ); im(1, 3) = det * ( m(0, 0) * A2312 - m(0, 2) * A0312 + m(0, 3) * A0212 ); - //im(2, 0) = det * ( m(1, 0) * A1323 - m(1, 1) * A0323 + m(1, 3) * A0123 ); - //im(2, 1) = det * - ( m(0, 0) * A1323 - m(0, 1) * A0323 + m(0, 3) * A0123 ); - //im(2, 2) = det * ( m(0, 0) * A1313 - m(0, 1) * A0313 + m(0, 3) * A0113 ); + // im(2, 0) = det * ( m(1, 0) * A1323 - m(1, 1) * A0323 + m(1, 3) * A0123 ); + // im(2, 1) = det * - ( m(0, 0) * A1323 - m(0, 1) * A0323 + m(0, 3) * A0123 ); + // im(2, 2) = det * ( m(0, 0) * A1313 - m(0, 1) * A0313 + m(0, 3) * A0113 ); im(2, 3) = det * - ( m(0, 0) * A1312 - m(0, 1) * A0312 + m(0, 3) * A0112 ); - //im(3, 0) = det * - ( m(1, 0) * A1223 - m(1, 1) * A0223 + m(1, 2) * A0123 ); - //im(3, 1) = det * ( m(0, 0) * A1223 - m(0, 1) * A0223 + m(0, 2) * A0123 ); - //im(3, 2) = det * - ( m(0, 0) * A1213 - m(0, 1) * A0213 + m(0, 2) * A0113 ); + // im(3, 0) = det * - ( m(1, 0) * A1223 - m(1, 1) * A0223 + m(1, 2) * A0123 ); + // im(3, 1) = det * ( m(0, 0) * A1223 - m(0, 1) * A0223 + m(0, 2) * A0123 ); + // im(3, 2) = det * - ( m(0, 0) * A1213 - m(0, 1) * A0213 + m(0, 2) * A0113 ); im(3, 3) = det * ( m(0, 0) * A1212 - m(0, 1) * A0212 + m(0, 2) * A0112 ); return true; } -template -Eigen::Matrix -vector_to_col_vec(const typename GeomTraits::Vector_3& v) +template +typename GarlandHeckbert_matrix_types::Col_3 +vector_to_col_3(const typename GeomTraits::Vector_3& v) { - Eigen::Matrix col {v.x(), v.y(), v.z()}; + typename GarlandHeckbert_matrix_types::Col_3 col { v.x(), v.y(), v.z() }; return col; } -// convenience alias declarations to make function return and argument types more readable -template -using Mat_4 = Eigen::Matrix; - -template -using Col_4 = Eigen::Matrix; - -template +template typename GeomTraits::Vector_3 construct_unit_normal_from_face(typename boost::graph_traits::face_descriptor f, const TriangleMesh& tmesh, - const VertexPointMap point_map, + const VertexPointMap vpm, const GeomTraits& gt) { - typedef typename boost::property_traits::reference point_reference; typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + typedef typename boost::property_traits::reference Point_ref; - auto unit_normal = gt.construct_unit_normal_3_object(); + typedef typename GeomTraits::FT FT; + typedef typename GeomTraits::Vector_3 Vector_3; + + auto cross_product = gt.construct_cross_product_vector_3_object(); + auto squared_length = gt.compute_squared_length_3_object(); const halfedge_descriptor h = halfedge(f, tmesh); - 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)); + const Point_ref p = get(vpm, target(h, tmesh)); + const Point_ref q = get(vpm, target(next(h, tmesh), tmesh)); + const Point_ref r = get(vpm, source(h, tmesh)); - return unit_normal(p, q, r); -} + Vector_3 normal = cross_product(q - p, r - p); -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; - - 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 = 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()); + const FT norm = sqrt(squared_length(normal)); + if(!is_zero(norm)) + normal = normal / norm; 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) +template +typename GeomTraits::Vector_3 +construct_edge_normal(typename boost::graph_traits::halfedge_descriptor h, + const TriangleMesh& tmesh, + const VertexPointMap vpm, + const GeomTraits& gt) { - typedef typename GeomTraits::FT FT; - - auto dot_product = gt.compute_scalar_product_3_object(); - 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)); - - // row vector given by d appended to the normal - 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(typename boost::graph_traits::face_descriptor f, - const TriangleMesh& mesh, - const VertexPointMap point_map, - const GeomTraits& gt) -{ - auto normal = construct_unit_normal_from_face(f, mesh, point_map, gt); - - // 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(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 boost::property_traits::reference Point_ref; - 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 - 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) -{ typedef typename GeomTraits::FT FT; typedef typename GeomTraits::Vector_3 Vector_3; - typedef Eigen::Matrix Mat_4; + auto vector = gt.construct_vector_3_object(); + auto cross_product = gt.construct_cross_product_vector_3_object(); 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); + const Point_ref p = get(vpm, target(h, tmesh)); + const Point_ref q = get(vpm, target(next(h, tmesh), tmesh)); + const Point_ref r = get(vpm, source(h, tmesh)); + const Vector_3 face_normal = cross_product(q - p, r - p); - // Eigen column vector of length 3 - const Eigen::Matrix mean_n_col { mean_normal.x(), mean_normal.y(), mean_normal.z() }; + const Vector_3 edge_vector = vector(r, p); + Vector_3 normal = cross_product(edge_vector, face_normal); - // start by setting values along the diagonal - Mat_4 mat = face_nv * Mat_4::Identity(); + const FT norm = sqrt(squared_length(normal)); + if(!is_zero(norm)) + normal = normal / norm; - // add outer product of the mean normal with itself - // to the upper left 3x3 block - mat.block(0, 0, 3, 3) += mean_n_col * mean_n_col.transpose(); - - // set the first 3 values of the last row and the first - // 3 values of the last column - // the negative sign comes from the fact that in the paper, - // 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() }; - - mat.col(3).head(3) = b; - mat.row(3).head(3) = b.transpose(); - - // 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; - - return mat; + return normal; } -template +template std::array -vectors_from_face(const VertexPointMap point_map, +vectors_from_face(typename boost::graph_traits::face_descriptor f, const TriangleMesh& tmesh, - typename boost::graph_traits::face_descriptor f, + const VertexPointMap vpm, 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(); + auto vector = gt.construct_vector_3_object(); const halfedge_descriptor h = halfedge(f, tmesh); // get all points and turn them into location vectors so we can use cross product on them - 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)); + const Point_reference p = get(vpm, source(h, tmesh)); + const Point_reference q = get(vpm, target(h, tmesh)); + const Point_reference r = get(vpm, target(next(h, tmesh), tmesh)); - std::array arr { construct_vector(ORIGIN, p), - construct_vector(ORIGIN, q), - construct_vector(ORIGIN, r) }; + std::array arr { vector(ORIGIN, p), + vector(ORIGIN, q), + vector(ORIGIN, r) }; return arr; } -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(); - - 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]; - - const Vector_3 ab = cross_product(a, b); - const Vector_3 bc = cross_product(b, c); - const Vector_3 ca = cross_product(c, a); - - const Vector_3 sum_of_cross_product = sum_vectors(sum_vectors(ab, bc), ca); - 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; - - // calculate the outer product of row^t*row - return row.transpose() * row; -} - -template -Eigen::Matrix +template +typename GarlandHeckbert_matrix_types::Mat_3 skew_sym_mat_cross_product(const typename GeomTraits::Vector_3& v) { - Eigen::Matrix mat; + typename GarlandHeckbert_matrix_types::Mat_3 mat; mat << 0, - v.z(), v.y(), v.z(), 0, - v.x(), @@ -316,114 +193,41 @@ skew_sym_mat_cross_product(const typename GeomTraits::Vector_3& v) return mat; } -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) +//////////////////////////////////////////////////////////////////////////////////////////////////// +/// OPTIMAL POINT +//////////////////////////////////////////////////////////////////////////////////////////////////// + +template +typename GarlandHeckbert_matrix_types::Col_4 +construct_optimal_point_invertible(const typename GarlandHeckbert_matrix_types::Mat_4& quadric) { - typedef typename GeomTraits::FT FT; - typedef typename GeomTraits::Vector_3 Vector_3; + typedef typename GarlandHeckbert_matrix_types::Mat_4 Mat_4; + typedef typename GarlandHeckbert_matrix_types::Col_4 Col_4; - 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); - - 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); - const Vector_3 bc = cross_product(b, c); - const Vector_3 ca = cross_product(c, a); - - // subtracting vectors using GeomTraits - const Vector_3 a_minus_b = sum_vectors(a, -b); - const Vector_3 b_minus_c = sum_vectors(b, -c); - const Vector_3 c_minus_a = sum_vectors(c, -a); - - const Mat_3 cp_ab = skew_sym_mat_cross_product(a_minus_b); - const Mat_3 cp_bc = skew_sym_mat_cross_product(b_minus_c); - const Mat_3 cp_ca = skew_sym_mat_cross_product(c_minus_a); - - 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() }; - - 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 * 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 * square(var) * vector_to_col_vec(sum_vectors(sum_vectors(a, b), c)); - - 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(3, 3) = res_c; - - return ret; -} - -template -Col_4 construct_optimal_point_invertible(const Mat_4& quadric) -{ - Mat_4 x; + Mat_4 x; x << quadric.block(0, 0, 3, 4), 0, 0, 0, 1; - Col_4 opt_pt; - + Col_4 opt_pt; opt_pt = x.inverse().col(3); // == X.inverse() * (0 0 0 1) + return opt_pt; } template -Col_4 construct_optimal_point_singular(const Mat_4& quadric, - const Col_4& p0, - const Col_4& p1) +typename GarlandHeckbert_matrix_types::Col_4 +construct_optimal_point_singular(const typename GarlandHeckbert_matrix_types::Mat_4& quadric, + const typename GarlandHeckbert_matrix_types::Col_4& p0, + const typename GarlandHeckbert_matrix_types::Col_4& p1) { typedef typename GeomTraits::FT FT; + typedef typename GarlandHeckbert_matrix_types::Mat_4 Mat_4; + typedef typename GarlandHeckbert_matrix_types::Col_4 Col_4; - // in this case, the matrix mat may no be invertible, - // so we save the result to check - Mat_4 mat; + // In this case, the matrix mat may no be invertible, so we save the result to check + Mat_4 mat; mat << quadric.block(0, 0, 3, 4), 0, 0, 0, 1; - Mat_4 inverse; + Mat_4 inverse; bool invertible = invert_matrix_4(mat, inverse); if(invertible) @@ -432,9 +236,9 @@ Col_4 construct_optimal_point_singular(const Mat_4& quad } else { - Col_4 opt_pt; + Col_4 opt_pt; - const Col_4 p1mp0 = p1 - p0; + const Col_4 p1mp0 = p1 - p0; const FT a = (p1mp0.transpose() * quadric * p1mp0)(0, 0); const FT b = 2 * (p0.transpose() * quadric * p1mp0)(0, 0); @@ -449,7 +253,7 @@ Col_4 construct_optimal_point_singular(const Mat_4& quad } else { - FT ext_t = -b / (2 * a); + FT ext_t = - b / (2 * a); if(ext_t < 0 || ext_t > 1 || a < 0) { // one of endpoints @@ -472,35 +276,288 @@ Col_4 construct_optimal_point_singular(const Mat_4& quad } } -template + +//////////////////////////////////////////////////////////////////////////////////////////////////// +/// CLASSIC PLANE +//////////////////////////////////////////////////////////////////////////////////////////////////// + +template +typename GarlandHeckbert_matrix_types::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; + + typedef typename GarlandHeckbert_matrix_types::Row_4 Row_4; + + auto dot_product = gt.compute_scalar_product_3_object(); + auto vector = gt.construct_vector_3_object(); + + // negative dot product between the normal and the position vector + const FT d = - dot_product(normal, vector(ORIGIN, point)); + + // row vector given by d appended to the normal + Row_4 row { normal.x(), normal.y(), normal.z(), d }; + + // outer product + return row.transpose() * row; +} + +template +typename GarlandHeckbert_matrix_types::Mat_4 +construct_classic_plane_quadric_from_edge(typename boost::graph_traits::halfedge_descriptor he, + const TriangleMesh& mesh, + const VertexPointMap vpm, + const GeomTraits& gt) +{ + 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, vpm, 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(vpm, source(he, mesh)), gt); +} + +template +typename GarlandHeckbert_matrix_types::Mat_4 +construct_classic_plane_quadric_from_face(typename boost::graph_traits::face_descriptor f, + const TriangleMesh& mesh, + const VertexPointMap vpm, + const GeomTraits& gt) +{ + auto normal = construct_unit_normal_from_face(f, mesh, vpm, gt); + + // get any point of the face + const auto p = get(vpm, source(halfedge(f, mesh), mesh)); + + return construct_classic_plane_quadric_from_normal(normal, p, gt); +} + +//////////////////////////////////////////////////////////////////////////////////////////////////// +/// CLASSIC TRIANGLE +//////////////////////////////////////////////////////////////////////////////////////////////////// + +template +typename GarlandHeckbert_matrix_types::Mat_4 +construct_classic_triangle_quadric_from_face(typename boost::graph_traits::face_descriptor f, + const TriangleMesh& tmesh, + const VertexPointMap vpm, + const GeomTraits& gt) +{ + typedef typename GeomTraits::FT FT; + typedef typename GeomTraits::Vector_3 Vector_3; + + typedef typename GarlandHeckbert_matrix_types::Row_4 Row_4; + + 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(); + + auto vectors = vectors_from_face(f, tmesh, vpm, gt); + + const Vector_3& a = vectors[0]; + const Vector_3& b = vectors[1]; + const Vector_3& c = vectors[2]; + + const Vector_3 ab = cross_product(a, b); + const Vector_3 bc = cross_product(b, c); + const Vector_3 ca = cross_product(c, a); + + const Vector_3 sum_of_cross_product = sum_vectors(sum_vectors(ab, bc), ca); + const FT scalar_triple_product = dot_product(ab, c); + + Row_4 row; + 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; +} + +//////////////////////////////////////////////////////////////////////////////////////////////////// +/// PROB PLANE +//////////////////////////////////////////////////////////////////////////////////////////////////// + +template +typename GarlandHeckbert_matrix_types::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 typename GarlandHeckbert_matrix_types::Col_3 Col_3; + typedef typename GarlandHeckbert_matrix_types::Mat_4 Mat_4; + + auto vector = gt.construct_vector_3_object(); + auto dot_product = gt.compute_scalar_product_3_object(); + auto squared_length = gt.compute_squared_length_3_object(); + + const Vector_3 mean_vec = vector(ORIGIN, point); + const FT dot_mnmv = dot_product(mean_normal, mean_vec); + + // Eigen column vector of length 3 + Col_3 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(); + + // add outer product of the mean normal with itself + // to the upper left 3x3 block + mat.block(0, 0, 3, 3) += mean_n_col * mean_n_col.transpose(); + + // set the first 3 values of the last row and the first + // 3 values of the last column + // the negative sign comes from the fact that in the paper, + // the b column and row appear with a negative sign + const auto b1 = -(dot_mnmv * mean_normal + face_nv * mean_vec); + + const Col_3 b { b1.x(), b1.y(), b1.z() }; + + mat.col(3).head(3) = b; + mat.row(3).head(3) = b.transpose(); + + // 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) = square(dot_mnmv) + + face_nv * squared_length(mean_vec) + + face_mv * squared_length(mean_normal) + + 3 * face_nv * face_mv; + + return mat; +} + +//////////////////////////////////////////////////////////////////////////////////////////////////// +/// PROB TRIANGLE +//////////////////////////////////////////////////////////////////////////////////////////////////// + +template +typename GarlandHeckbert_matrix_types::Mat_4 +construct_prob_triangle_quadric_from_face(typename boost::graph_traits::face_descriptor f, + typename GeomTraits::FT var, + const TriangleMesh& tmesh, + const VertexPointMap vpm, + const GeomTraits& gt) +{ + typedef typename GeomTraits::FT FT; + typedef typename GeomTraits::Vector_3 Vector_3; + + typedef typename GarlandHeckbert_matrix_types::Mat_3 Mat_3; + typedef typename GarlandHeckbert_matrix_types::Col_3 Col_3; + typedef typename GarlandHeckbert_matrix_types::Mat_4 Mat_4; + + 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(f, tmesh, vpm, gt); + + 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); + const Vector_3 bc = cross_product(b, c); + const Vector_3 ca = cross_product(c, a); + + // subtracting vectors using GeomTraits + const Vector_3 a_minus_b = sum_vectors(a, -b); + const Vector_3 b_minus_c = sum_vectors(b, -c); + const Vector_3 c_minus_a = sum_vectors(c, -a); + + const Mat_3 cp_ab = skew_sym_mat_cross_product(a_minus_b); + const Mat_3 cp_bc = skew_sym_mat_cross_product(b_minus_c); + const Mat_3 cp_ca = skew_sym_mat_cross_product(c_minus_a); + + const Vector_3 sum_of_cross_product = sum_vectors(sum_vectors(ab, bc), ca); + + const Col_3 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 * 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 Col_3 res_b = det * sum_cp_col - var * (vector_to_col_3(cross_product(a_minus_b, ab)) + + vector_to_col_3(cross_product(b_minus_c, bc)) + + vector_to_col_3(cross_product(c_minus_a, ca))) + + 2 * square(var) * vector_to_col_3(sum_vectors(sum_vectors(a, b), c)); + + 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(3, 3) = res_c; + + return ret; +} + +//////////////////////////////////////////////////////////////////////////////////////////////////// +/// PROB VARIANCE +//////////////////////////////////////////////////////////////////////////////////////////////////// + +template std::pair estimate_variances(const TriangleMesh& mesh, const GeomTraits& gt, - typename GeomTraits::FT variance, - typename GeomTraits::FT p_factor) + const typename GeomTraits::FT variance, + const 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 GeomTraits::Vector_3 Vector_3; + + auto construct_vector = gt.construct_vector_3_object(); + auto squared_length = gt.compute_squared_length_3_object(); FT average_edge_length = 0; - auto construct_vector = gt.construct_vector_3_object(); - + std::size_t ne = 0; 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& p1 = mesh.point(v1); // @fixme Surface_mesh API const Point_3& p2 = mesh.point(v2); const Vector_3 vec = construct_vector(p1, p2); - average_edge_length += sqrt(vec.squared_length()); + average_edge_length += sqrt(squared_length(vec)); + + ++ne; // edges(mesh).size() can be costly, might as well increment now } - average_edge_length = average_edge_length / mesh.number_of_edges(); + average_edge_length = average_edge_length / ne; + const FT n2 = variance * average_edge_length; const FT p2 = p_factor * variance * average_edge_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 dfe6156d28c..e4f40ec503a 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 @@ -28,149 +28,138 @@ namespace CGAL { namespace Surface_mesh_simplification { namespace internal { -template -Mat_4 combine_matrices(const Mat_4& a, const Mat_4& b) +template +Mat combine_matrices(const Mat& a, const Mat& b) { return a + b; } -template -class GarlandHeckbert_cost_base +template +struct GarlandHeckbert_matrix_types { + typedef typename GeomTraits::FT FT; + + typedef Eigen::Matrix Mat_3; + typedef Eigen::Matrix Col_3; + typedef Eigen::Matrix Mat_4; + typedef Eigen::Matrix Col_4; + typedef Eigen::Matrix Row_4; +}; + +// Storage is initialized by the most-derived class (e.g. GarlandHeckbert_plane_policies) +template +struct GarlandHeckbert_quadrics_storage +{ + typedef typename GarlandHeckbert_matrix_types::Mat_4 Mat_4; + typedef typename GarlandHeckbert_matrix_types::Col_4 Col_4; + + typedef Mat_4 Cost_matrix; + typedef CGAL::dynamic_vertex_property_t Cost_property; + typedef typename boost::property_map::type Vertex_cost_map; + + typedef QuadricCalculator Quadric_calculator; + +protected: + Vertex_cost_map m_cost_matrices; + Quadric_calculator m_quadric_calculator; + +public: + GarlandHeckbert_quadrics_storage() { CGAL_assertion(false); } // only for compilation + GarlandHeckbert_quadrics_storage(TriangleMesh& tmesh) + { + m_cost_matrices = get(Cost_property(), tmesh); + } + + GarlandHeckbert_quadrics_storage(TriangleMesh& tmesh, + const Quadric_calculator& quadric_calculator) + : m_quadric_calculator(quadric_calculator) + { + m_cost_matrices = get(Cost_property(), tmesh); + } +}; + +template +class GarlandHeckbert_cost_base + // Virtual base because the storage is common to both the cost and the placement + : virtual public GarlandHeckbert_quadrics_storage +{ + typedef QuadricCalculator Quadric_calculator; + typedef GarlandHeckbert_quadrics_storage< + Quadric_calculator, TriangleMesh, GeomTraits> 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 VertexCostMap Vertex_cost_map; + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + typedef typename boost::graph_traits::face_descriptor face_descriptor; + typedef typename GeomTraits::FT FT; - - typedef Eigen::Matrix Mat_4; - typedef Eigen::Matrix Col_4; - typedef typename GeomTraits::Point_3 Point_3; typedef typename GeomTraits::Vector_3 Vector_3; -private: - Vertex_cost_map m_cost_matrices; + typedef typename Base::Mat_4 Mat_4; + typedef typename Base::Col_4 Col_4; +private: FT discontinuity_multiplier; public: - GarlandHeckbert_cost_base() - : m_cost_matrices() + GarlandHeckbert_cost_base(const FT dm = 100) + : discontinuity_multiplier(dm) { } - GarlandHeckbert_cost_base(FT dm) - : m_cost_matrices(), discontinuity_multiplier(dm) - { } + decltype(auto) vcm() const { return this->m_cost_matrices; } + const Quadric_calculator& quadric_calculator() const { return this->m_quadric_calculator; } - GarlandHeckbert_cost_base(Vertex_cost_map vcm, FT dm) - : m_cost_matrices(vcm), discontinuity_multiplier(dm) - { } - -protected: - void init_vcm(const Vertex_cost_map vcm) +public: + static Col_4 point_to_homogenous_column(const Point_3& point) { - m_cost_matrices = vcm; + return Col_4 { point.x(), point.y(), point.z(), FT(1) }; } - 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) + static bool is_discontinuity_edge(const halfedge_descriptor h, + const TriangleMesh& tmesh) { return is_border_edge(h, tmesh); } +public: + template + Mat_4 construct_quadric(const halfedge_descriptor he, + const TriangleMesh& tmesh, + const VertexPointMap vpm, + const GeomTraits& gt) const + { + return quadric_calculator().construct_quadric_from_edge(he, tmesh, vpm, gt); + } + + template + Mat_4 construct_quadric(const face_descriptor f, + const TriangleMesh& tmesh, + const VertexPointMap vpm, + const GeomTraits& gt) const + { + return quadric_calculator().construct_quadric_from_face(f, tmesh, vpm, gt); + } + public: // initialize all quadrics - template - void initialize(const TM& tmesh, - const VPM vpm, + template + void initialize(const TriangleMesh& tmesh, + const VertexPointMap vpm, const GeomTraits& gt) const { - 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); + put(vcm(), v, zero_mat); for(face_descriptor f : faces(tmesh)) { - if(f == GraphTraits::null_face()) + if(f == boost::graph_traits::null_face()) continue; const halfedge_descriptor h = halfedge(f, tmesh); @@ -183,7 +172,7 @@ public: const vertex_descriptor vs = source(shd, tmesh); const vertex_descriptor vt = target(shd, tmesh); - put(m_cost_matrices, vs, combine_matrices(get(m_cost_matrices, vs), quadric)); + put(vcm(), vs, combine_matrices(get(vcm(), vs), quadric)); if(!is_discontinuity_edge(shd, tmesh)) continue; @@ -191,12 +180,21 @@ public: 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)); + put(vcm(), vs, combine_matrices(get(vcm(), vs), discontinuous_quadric)); + put(vcm(), vt, combine_matrices(get(vcm(), vt), discontinuous_quadric)); } } } + template + void update_after_collapse(const Profile& profile, + const VertexDescriptor new_v) const + { + put(vcm(), new_v, combine_matrices(get(vcm(), profile.v0()), + get(vcm(), profile.v1()))); + } + +public: template boost::optional operator()(const Profile& profile, @@ -205,88 +203,69 @@ public: typedef boost::optional Optional_FT; if(!placement) - { - // return empty value return boost::optional(); - } - CGAL_precondition(!get(m_cost_matrices, profile.v0()).isZero(0)); - CGAL_precondition(!get(m_cost_matrices, profile.v1()).isZero(0)); + CGAL_precondition(!get(vcm(), profile.v0()).isZero(0)); + CGAL_precondition(!get(vcm(), profile.v1()).isZero(0)); - const Mat_4 combined_matrix = combine_matrices(get(m_cost_matrices, profile.v0()), - get(m_cost_matrices, profile.v1())); + const Mat_4 combined_matrix = combine_matrices(get(vcm(), profile.v0()), + get(vcm(), profile.v1())); const Col_4 pt = point_to_homogenous_column(*placement); const Optional_FT cost = (pt.transpose() * combined_matrix * pt)(0, 0); return cost; } - - template - void update_after_collapse(const Profile& profile, - 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()))); - } }; -template +template class GarlandHeckbert_placement_base + // Virtual base because the storage is common to both the cost and the placement + : virtual public GarlandHeckbert_quadrics_storage { -public: - // define type required by the Get_cost concept - typedef VertexCostMap Vertex_cost_map; + typedef GarlandHeckbert_quadrics_storage< + QuadricCalculator, TriangleMesh, GeomTraits> Base; - // matrix and column vector types +public: typedef typename GeomTraits::FT FT; typedef typename GeomTraits::Point_3 Point_3; - typedef Eigen::Matrix Mat_4; - typedef Eigen::Matrix Col_4; -private: - Vertex_cost_map m_cost_matrices; + typedef QuadricCalculator Quadric_calculator; + + typedef typename GarlandHeckbert_matrix_types::Mat_4 Mat_4; + typedef typename GarlandHeckbert_matrix_types::Col_4 Col_4; public: GarlandHeckbert_placement_base() { } - GarlandHeckbert_placement_base(Vertex_cost_map cost_matrices) - : m_cost_matrices(cost_matrices) - { } + + decltype(auto) vcm() const { return this->m_cost_matrices; } + const Quadric_calculator& quadric_calculator() const { return this->m_quadric_calculator; } protected: - void init_vcm(const Vertex_cost_map& vcm) + static Col_4 point_to_homogenous_column(const Point_3& point) { - m_cost_matrices = vcm; + return Col_4 { point.x(), point.y(), point.z(), FT(1) }; } - // use CRTP to call the quadric implementation 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 - { - return Col_4(point.x(), point.y(), point.z(), FT(1)); + return quadric_calculator().construct_optimal_point(mat, p0, p1); } 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)); + CGAL_precondition(!get(vcm(), profile.v0()).isZero(0)); + CGAL_precondition(!get(vcm(), 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 Mat_4 combined_matrix = combine_matrices(get(vcm(), profile.v0()), + get(vcm(), 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); + const Col_4 opt = construct_optimum(combined_matrix, p0, p1); boost::optional pt = typename Profile::Point(opt(0) / opt(3), opt(1) / opt(3),