songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Refactor + improvements for Garland Heckbert methods

This commit is contained in:
Mael Rouxel-Labbé 2021-12-22 12:03:11 +01:00
parent a1b00f9206
commit b342e11967
6 changed files with 848 additions and 845 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

153
Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_policies.h
View File

@ -16,95 +16,22 @@
#include <CGAL/license/Surface_mesh_simplification.h>
#include <CGAL/Surface_mesh_simplification/internal/Common.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/GarlandHeckbert_functions.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/GarlandHeckbert_policy_base.h>
#include <Eigen/Dense>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/GarlandHeckbert_functions.h>
namespace CGAL {
namespace Surface_mesh_simplification {
template<typename TriangleMesh, typename GeomTraits>
class GarlandHeckbert_policies :
public internal::GarlandHeckbert_placement_base<
typename boost::property_map<
TriangleMesh,
CGAL::dynamic_vertex_property_t<Eigen::Matrix<typename GeomTraits::FT, 4, 4, Eigen::DontAlign> >
>::type,
GeomTraits,
GarlandHeckbert_policies<TriangleMesh, GeomTraits>
>,
public internal::GarlandHeckbert_cost_base<
typename boost::property_map<
TriangleMesh,
CGAL::dynamic_vertex_property_t<Eigen::Matrix<typename GeomTraits::FT, 4, 4, Eigen::DontAlign> >
>::type,
GeomTraits,
GarlandHeckbert_policies<TriangleMesh, GeomTraits> >
template <typename TriangleMesh, typename GeomTraits>
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<FT, 4, 4, Eigen::DontAlign> GH_matrix;
typedef CGAL::dynamic_vertex_property_t<GH_matrix> Cost_property;
typedef typename boost::property_map<TriangleMesh, Cost_property>::type Vertex_cost_map;
typedef internal::GarlandHeckbert_placement_base<
Vertex_cost_map, GeomTraits, GarlandHeckbert_policies<TriangleMesh, GeomTraits>
> Placement_base;
typedef internal::GarlandHeckbert_cost_base<
Vertex_cost_map, GeomTraits, GarlandHeckbert_policies<TriangleMesh, GeomTraits>
> Cost_base;
// both types are the same, this is so we avoid casting back to the base class in 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<GeomTraits>::Mat_4 Mat_4;
typedef typename internal::GarlandHeckbert_matrix_types<GeomTraits>::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<typename VertexPointMap>
Mat_4 construct_quadric_from_face(typename boost::graph_traits<TriangleMesh>::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<typename VertexPointMap>
template <typename VertexPointMap>
Mat_4 construct_quadric_from_edge(typename boost::graph_traits<TriangleMesh>::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 <typename VertexPointMap>
Mat_4 construct_quadric_from_face(typename boost::graph_traits<TriangleMesh>::face_descriptor f,
const TriangleMesh& tmesh,
const VertexPointMap point_map,
const GeomTraits& gt) const
{
return internal::construct_optimal_point_singular<GeomTraits>(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<FT, 1, 4> 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<GeomTraits>(quadric, p0, p1);
}
};
template<typename TriangleMesh, typename GeomTraits>
class GarlandHeckbert_policies
: public internal::GarlandHeckbert_placement_base<
Plane_quadric_calculator<TriangleMesh, GeomTraits>, TriangleMesh, GeomTraits>,
public internal::GarlandHeckbert_cost_base<
Plane_quadric_calculator<TriangleMesh, GeomTraits>, TriangleMesh, GeomTraits>
{
public:
typedef Plane_quadric_calculator<TriangleMesh, GeomTraits> 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<TriangleMesh, GeomTraits> 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

212
Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_probabilistic_policies.h
View File

@ -19,161 +19,165 @@
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/GarlandHeckbert_policy_base.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/GarlandHeckbert_functions.h>
#include <Eigen/Dense>
#include <CGAL/Default.h>
#include <CGAL/property_map.h>
#include <tuple>
#include <utility>
namespace CGAL {
namespace Surface_mesh_simplification {
// Implements probabilistic plane quadrics,
// optionally takes a face variance map giving a per-face variance
template<typename TriangleMesh,
typename GeomTraits,
typename FaceVarianceMap = Constant_property_map<
typename boost::graph_traits<TriangleMesh>::face_descriptor,
std::pair<typename GeomTraits::FT, typename GeomTraits::FT> > >
class GarlandHeckbert_probabilistic_policies
: public internal::GarlandHeckbert_placement_base<
typename boost::property_map<
TriangleMesh,
CGAL::dynamic_vertex_property_t<Eigen::Matrix<typename GeomTraits::FT, 4, 4, Eigen::DontAlign> >
>::type,
GeomTraits,
GarlandHeckbert_probabilistic_policies<TriangleMesh, GeomTraits, FaceVarianceMap> >,
public internal::GarlandHeckbert_cost_base<
typename boost::property_map<
TriangleMesh,
CGAL::dynamic_vertex_property_t<Eigen::Matrix<typename GeomTraits::FT, 4, 4, Eigen::DontAlign>>
>::type,
GeomTraits,
GarlandHeckbert_probabilistic_policies<TriangleMesh, GeomTraits, FaceVarianceMap> >
template <typename TriangleMesh, typename GeomTraits, typename FaceVarianceMap>
class Probabilistic_plane_quadric_calculator
{
typedef typename GeomTraits::FT FT;
typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::graph_traits<TriangleMesh>::face_descriptor face_descriptor;
typedef typename boost::property_traits<FaceVarianceMap>::value_type Face_variance;
typedef typename GeomTraits::FT FT;
public:
typedef typename Eigen::Matrix<FT, 4, 4, Eigen::DontAlign> GH_matrix;
typedef CGAL::dynamic_vertex_property_t<GH_matrix> Cost_property;
typedef typename boost::property_map<TriangleMesh, Cost_property>::type Vertex_cost_map;
typedef Constant_property_map<face_descriptor, std::pair<FT, FT> > Default_FVM;
typedef typename Default::Get<FaceVarianceMap, Default_FVM>::type Face_variance_map;
typedef internal::GarlandHeckbert_placement_base<
Vertex_cost_map, GeomTraits, GarlandHeckbert_probabilistic_policies<TriangleMesh, GeomTraits>
> Placement_base;
typedef internal::GarlandHeckbert_cost_base<
Vertex_cost_map, GeomTraits, GarlandHeckbert_probabilistic_policies<TriangleMesh, GeomTraits>
> 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<GeomTraits>::Mat_4 Mat_4;
typedef typename internal::GarlandHeckbert_matrix_types<GeomTraits>::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 <typename FVM>
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<std::is_same<Face_variance_map, Default_FVM> >::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<typename VPM>
Mat_4 construct_quadric_from_face(face_descriptor f,
template <typename VertexPointMap>
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<typename VPM>
Mat_4 construct_quadric_from_edge(halfedge_descriptor he,
template <typename VertexPointMap>
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<GeomTraits>(quadric);
}
};
// Implements probabilistic plane quadrics,
// optionally takes a face variance map giving a per-face variance
template<typename TriangleMesh,
typename GeomTraits,
typename FaceVarianceMap = CGAL::Default>
class GarlandHeckbert_probabilistic_policies
: public internal::GarlandHeckbert_placement_base<
Probabilistic_plane_quadric_calculator<TriangleMesh, GeomTraits, FaceVarianceMap>,
TriangleMesh, GeomTraits>,
public internal::GarlandHeckbert_cost_base<
Probabilistic_plane_quadric_calculator<TriangleMesh, GeomTraits, FaceVarianceMap>,
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<TriangleMesh, GeomTraits> 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 <typename FVM>
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

220
Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_probabilistic_tri_policies.h
View File

@ -17,167 +17,155 @@
#include <CGAL/Surface_mesh_simplification/internal/Common.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/GarlandHeckbert_policy_base.h>
#include <Eigen/Dense>
#include <CGAL/Default.h>
#include <CGAL/property_map.h>
#include <CGAL/boost/graph/Named_function_parameters.h>
#include <utility>
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<typename TriangleMesh,
typename GeomTraits,
typename FaceVarianceMap =
Constant_property_map<typename boost::graph_traits<TriangleMesh>::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<Eigen::Matrix<typename GeomTraits::FT, 4, 4, Eigen::DontAlign>>
>::type,
GeomTraits,
GarlandHeckbert_probabilistic_tri_policies<TriangleMesh, GeomTraits, FaceVarianceMap>
>,
public internal::GarlandHeckbert_cost_base<
typename boost::property_map<
TriangleMesh,
CGAL::dynamic_vertex_property_t<Eigen::Matrix<typename GeomTraits::FT, 4, 4, Eigen::DontAlign>>
>::type,
GeomTraits,
GarlandHeckbert_probabilistic_tri_policies<TriangleMesh, GeomTraits, FaceVarianceMap> >
template <typename TriangleMesh, typename GeomTraits, typename FaceVarianceMap>
class Probabilistic_triangle_quadric_calculator
{
typedef typename boost::property_traits<FaceVarianceMap>::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<FT, 3, 3, Eigen::DontAlign> Mat_3;
typedef typename Eigen::Matrix<FT, 4, 4, Eigen::DontAlign> GH_matrix;
typedef CGAL::dynamic_vertex_property_t<GH_matrix> Cost_property;
typedef typename boost::property_map<TriangleMesh, Cost_property>::type Vertex_cost_map;
typedef typename boost::graph_traits<TriangleMesh>::face_descriptor face_descriptor;
typedef internal::GarlandHeckbert_placement_base<
Vertex_cost_map, GeomTraits, GarlandHeckbert_probabilistic_tri_policies<TriangleMesh, GeomTraits>
> Placement_base;
typedef Constant_property_map<face_descriptor, FT> Default_FVM;
typedef typename Default::Get<FaceVarianceMap, Default_FVM>::type Face_variance_map;
typedef internal::GarlandHeckbert_cost_base<
Vertex_cost_map, GeomTraits, GarlandHeckbert_probabilistic_tri_policies<TriangleMesh, GeomTraits>
> 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<GeomTraits>::Mat_4 Mat_4;
typedef typename internal::GarlandHeckbert_matrix_types<GeomTraits>::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 <typename FVM>
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<std::is_same<Face_variance_map, Default_FVM> >::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<typename VPM, typename TM>
Mat_4 construct_quadric_from_face(typename boost::graph_traits<TM>::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<typename VPM, typename TM>
Mat_4 construct_quadric_from_edge(typename boost::graph_traits<TM>::halfedge_descriptor he,
const TM& tmesh,
const VPM point_map,
template<typename VertexPointMap>
Mat_4 construct_quadric_from_edge(typename boost::graph_traits<TriangleMesh>::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<typename VertexPointMap>
Mat_4 construct_quadric_from_face(typename boost::graph_traits<TriangleMesh>::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<GeomTraits>(quadric);
}
};
// implements probabilistic triangle faces and optionally takes a face variance map
// analogously to probabilistic plane quadrics
template<typename TriangleMesh,
typename GeomTraits,
typename FaceVarianceMap = CGAL::Default>
class GarlandHeckbert_probabilistic_tri_policies
: public internal::GarlandHeckbert_placement_base<
Probabilistic_triangle_quadric_calculator<TriangleMesh, GeomTraits, FaceVarianceMap>,
TriangleMesh, GeomTraits>,
public internal::GarlandHeckbert_cost_base<
Probabilistic_triangle_quadric_calculator<TriangleMesh, GeomTraits, FaceVarianceMap>,
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<TriangleMesh, GeomTraits> 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 <typename FVM>
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

134
Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_triangle_policies.h
View File

@ -19,76 +19,29 @@
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/GarlandHeckbert_policy_base.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/GarlandHeckbert_functions.h>
#include <Eigen/Dense>
#include <iostream>
namespace CGAL {
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<typename TriangleMesh, typename GeomTraits>
class GarlandHeckbert_triangle_policies
: public internal::GarlandHeckbert_cost_base<
typename boost::property_map<
TriangleMesh,
CGAL::dynamic_vertex_property_t<Eigen::Matrix<typename GeomTraits::FT, 4, 4, Eigen::DontAlign>>
>::type,
GeomTraits,
GarlandHeckbert_triangle_policies<TriangleMesh, GeomTraits>
>,
public internal::GarlandHeckbert_placement_base<
typename boost::property_map<
TriangleMesh,
CGAL::dynamic_vertex_property_t<Eigen::Matrix<typename GeomTraits::FT, 4, 4, Eigen::DontAlign>>
>::type,
GeomTraits,
GarlandHeckbert_triangle_policies<TriangleMesh, GeomTraits> >
template <typename TriangleMesh, typename GeomTraits>
class Triangle_quadric_calculator
{
public:
typedef typename GeomTraits::FT FT;
typedef Eigen::Matrix<FT, 4, 4, Eigen::DontAlign> Mat_4;
typedef Eigen::Matrix<FT, 1, 4> Col_4;
typedef CGAL::dynamic_vertex_property_t<Mat_4> Cost_property;
typedef typename boost::property_map<TriangleMesh, Cost_property>::type Vertex_cost_map;
typedef internal::GarlandHeckbert_placement_base<
Vertex_cost_map, GeomTraits, GarlandHeckbert_triangle_policies<TriangleMesh, GeomTraits>
> Placement_base;
typedef internal::GarlandHeckbert_cost_base<
Vertex_cost_map, GeomTraits, GarlandHeckbert_triangle_policies<TriangleMesh, GeomTraits>
> 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<GeomTraits>::Mat_4 Mat_4;
typedef typename internal::GarlandHeckbert_matrix_types<GeomTraits>::Col_4 Col_4;
public:
GarlandHeckbert_triangle_policies(TriangleMesh& tmesh, FT dm = FT(100))
Triangle_quadric_calculator() { }
template <typename VertexPointMap>
Mat_4 construct_quadric_from_edge(typename boost::graph_traits<TriangleMesh>::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<typename VertexPointMap>
template <typename VertexPointMap>
Mat_4 construct_quadric_from_face(typename boost::graph_traits<TriangleMesh>::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<typename VertexPointMap>
Mat_4 construct_quadric_from_edge(typename boost::graph_traits<TriangleMesh>::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<GeomTraits>(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<typename TriangleMesh, typename GeomTraits>
class GarlandHeckbert_triangle_policies
: public internal::GarlandHeckbert_placement_base<
Triangle_quadric_calculator<TriangleMesh, GeomTraits>, TriangleMesh, GeomTraits>,
public internal::GarlandHeckbert_cost_base<
Triangle_quadric_calculator<TriangleMesh, GeomTraits>, TriangleMesh, GeomTraits>
{
public:
typedef Triangle_quadric_calculator<TriangleMesh, GeomTraits> 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<TriangleMesh, GeomTraits> 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

691
Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/GarlandHeckbert_functions.h
View File

@ -16,7 +16,7 @@
#include <CGAL/Surface_mesh_simplification/internal/Common.h>
#include <Eigen/Dense>
#include <CGAL/Origin.h>
#include <iostream>
@ -25,289 +25,166 @@ namespace Surface_mesh_simplification {
namespace internal {
// taken from https://stackoverflow.com/questions/1148309/inverting-a-4x4-matrix/
template<typename Matrix>
template <typename Matrix>
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<typename GeomTraits>
Eigen::Matrix<typename GeomTraits::FT, 3, 1>
vector_to_col_vec(const typename GeomTraits::Vector_3& v)
template <typename GeomTraits>
typename GarlandHeckbert_matrix_types<GeomTraits>::Col_3
vector_to_col_3(const typename GeomTraits::Vector_3& v)
{
Eigen::Matrix<typename GeomTraits::FT, 3, 1> col {v.x(), v.y(), v.z()};
typename GarlandHeckbert_matrix_types<GeomTraits>::Col_3 col { v.x(), v.y(), v.z() };
return col;
}
// convenience alias declarations to make function return and argument types more readable
template<typename GeomTraits>
using Mat_4 = Eigen::Matrix<typename GeomTraits::FT, 4, 4, Eigen::DontAlign>;
template<typename GeomTraits>
using Col_4 = Eigen::Matrix<typename GeomTraits::FT, 4, 1>;
template<typename VertexPointMap, typename TriangleMesh, typename GeomTraits>
template <typename TriangleMesh, typename VertexPointMap, typename GeomTraits>
typename GeomTraits::Vector_3
construct_unit_normal_from_face(typename boost::graph_traits<TriangleMesh>::face_descriptor f,
const TriangleMesh& tmesh,
const VertexPointMap point_map,
const VertexPointMap vpm,
const GeomTraits& gt)
{
typedef typename boost::property_traits<VertexPointMap>::reference point_reference;
typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::property_traits<VertexPointMap>::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 TriangleMesh, typename VertexPointMap, typename GeomTraits>
typename GeomTraits::Vector_3
construct_edge_normal(typename boost::graph_traits<TriangleMesh>::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<TriangleMesh>::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<typename GeomTraits>
Mat_4<GeomTraits>
construct_classic_plane_quadric_from_normal(const typename GeomTraits::Vector_3& normal,
const typename GeomTraits::Point_3& point,
const GeomTraits& gt)
template <typename TriangleMesh, typename VertexPointMap, typename GeomTraits>
typename GeomTraits::Vector_3
construct_edge_normal(typename boost::graph_traits<TriangleMesh>::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<FT, 1, 4> row { normal.x(), normal.y(), normal.z(), d };
// outer product
return row.transpose() * row;
}
template<typename TriangleMesh, typename VertexPointMap, typename GeomTraits>
Mat_4<GeomTraits>
construct_classic_plane_quadric_from_face(typename boost::graph_traits<TriangleMesh>::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<typename TriangleMesh, typename VertexPointMap, typename GeomTraits>
Mat_4<GeomTraits>
construct_classic_plane_quadric_from_edge(typename boost::graph_traits<TriangleMesh>::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<TriangleMesh>::vertex_descriptor vertex_descriptor;
typedef typename boost::property_traits<VertexPointMap>::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 <typename GeomTraits>
Mat_4<GeomTraits>
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<FT, 4, 4> 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<FT, 3, 1> 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<FT, 3, 1> 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<typename VertexPointMap, typename TriangleMesh, typename GeomTraits>
template <typename VertexPointMap, typename TriangleMesh, typename GeomTraits>
std::array<typename GeomTraits::Vector_3, 3>
vectors_from_face(const VertexPointMap point_map,
vectors_from_face(typename boost::graph_traits<TriangleMesh>::face_descriptor f,
const TriangleMesh& tmesh,
typename boost::graph_traits<TriangleMesh>::face_descriptor f,
const VertexPointMap vpm,
const GeomTraits& gt)
{
typedef typename boost::property_traits<VertexPointMap>::reference Point_reference;
typedef typename boost::graph_traits<TriangleMesh>::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<typename GeomTraits::Vector_3, 3> arr { construct_vector(ORIGIN, p),
construct_vector(ORIGIN, q),
construct_vector(ORIGIN, r) };
std::array<typename GeomTraits::Vector_3, 3> arr { vector(ORIGIN, p),
vector(ORIGIN, q),
vector(ORIGIN, r) };
return arr;
}
template<typename TriangleMesh, typename VertexPointMap, typename GeomTraits>
Mat_4<GeomTraits>
construct_classic_triangle_quadric_from_face(typename boost::graph_traits<TriangleMesh>::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<FT, 1, 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;
}
template<typename GeomTraits>
Eigen::Matrix<typename GeomTraits::FT, 3, 3>
template <typename GeomTraits>
typename GarlandHeckbert_matrix_types<GeomTraits>::Mat_3
skew_sym_mat_cross_product(const typename GeomTraits::Vector_3& v)
{
Eigen::Matrix<typename GeomTraits::FT, 3, 3> mat;
typename GarlandHeckbert_matrix_types<GeomTraits>::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<typename TriangleMesh, typename VertexPointMap, typename GeomTraits>
Mat_4<GeomTraits>
construct_prob_triangle_quadric_from_face(typename boost::graph_traits<TriangleMesh>::face_descriptor f,
typename GeomTraits::FT var,
const TriangleMesh& tmesh,
const VertexPointMap point_map,
const GeomTraits& gt)
////////////////////////////////////////////////////////////////////////////////////////////////////
/// OPTIMAL POINT
////////////////////////////////////////////////////////////////////////////////////////////////////
template <typename GeomTraits>
typename GarlandHeckbert_matrix_types<GeomTraits>::Col_4
construct_optimal_point_invertible(const typename GarlandHeckbert_matrix_types<GeomTraits>::Mat_4& quadric)
{
typedef typename GeomTraits::FT FT;
typedef typename GeomTraits::Vector_3 Vector_3;
typedef typename GarlandHeckbert_matrix_types<GeomTraits>::Mat_4 Mat_4;
typedef typename GarlandHeckbert_matrix_types<GeomTraits>::Col_4 Col_4;
typedef Eigen::Matrix<FT, 3, 3> Mat_3;
typedef Mat_4<GeomTraits> Mat_4;
auto construct_vector = gt.construct_vector_3_object();
auto cross_product = gt.construct_cross_product_vector_3_object();
auto sum_vectors = gt.construct_sum_of_vectors_3_object();
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<GeomTraits>(a_minus_b);
const Mat_3 cp_bc = skew_sym_mat_cross_product<GeomTraits>(b_minus_c);
const Mat_3 cp_ca = skew_sym_mat_cross_product<GeomTraits>(c_minus_a);
const Vector_3 sum_of_cross_product = sum_vectors(sum_vectors(ab, bc), ca);
const Eigen::Matrix<FT, 3, 1, Eigen::DontAlign>
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<FT, 3, 1> res_b =
det * sum_cp_col - var * (vector_to_col_vec<GeomTraits>(cross_product(a_minus_b, ab))
+ vector_to_col_vec<GeomTraits>(cross_product(b_minus_c, bc))
+ vector_to_col_vec<GeomTraits>(cross_product(c_minus_a, ca)))
+ 2 * square(var) * vector_to_col_vec<GeomTraits>(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<typename GeomTraits>
Col_4<GeomTraits> construct_optimal_point_invertible(const Mat_4<GeomTraits>& quadric)
{
Mat_4<GeomTraits> x;
Mat_4 x;
x << quadric.block(0, 0, 3, 4), 0, 0, 0, 1;
Col_4<GeomTraits> opt_pt;
Col_4 opt_pt;
opt_pt = x.inverse().col(3); // == X.inverse() * (0 0 0 1)
return opt_pt;
}
template <typename GeomTraits>
Col_4<GeomTraits> construct_optimal_point_singular(const Mat_4<GeomTraits>& quadric,
const Col_4<GeomTraits>& p0,
const Col_4<GeomTraits>& p1)
typename GarlandHeckbert_matrix_types<GeomTraits>::Col_4
construct_optimal_point_singular(const typename GarlandHeckbert_matrix_types<GeomTraits>::Mat_4& quadric,
const typename GarlandHeckbert_matrix_types<GeomTraits>::Col_4& p0,
const typename GarlandHeckbert_matrix_types<GeomTraits>::Col_4& p1)
{
typedef typename GeomTraits::FT FT;
typedef typename GarlandHeckbert_matrix_types<GeomTraits>::Mat_4 Mat_4;
typedef typename GarlandHeckbert_matrix_types<GeomTraits>::Col_4 Col_4;
// in this case, the matrix mat may no be invertible,
// so we save the result to check
Mat_4<GeomTraits> 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<GeomTraits> inverse;
Mat_4 inverse;
bool invertible = invert_matrix_4(mat, inverse);
if(invertible)
@ -432,9 +236,9 @@ Col_4<GeomTraits> construct_optimal_point_singular(const Mat_4<GeomTraits>& quad
}
else
{
Col_4<GeomTraits> opt_pt;
Col_4 opt_pt;
const Col_4<GeomTraits> 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<GeomTraits> construct_optimal_point_singular(const Mat_4<GeomTraits>& 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<GeomTraits> construct_optimal_point_singular(const Mat_4<GeomTraits>& quad
}
}
template<typename TriangleMesh, typename GeomTraits>
////////////////////////////////////////////////////////////////////////////////////////////////////
/// CLASSIC PLANE
////////////////////////////////////////////////////////////////////////////////////////////////////
template <typename GeomTraits>
typename GarlandHeckbert_matrix_types<GeomTraits>::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<GeomTraits>::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 TriangleMesh, typename VertexPointMap, typename GeomTraits>
typename GarlandHeckbert_matrix_types<GeomTraits>::Mat_4
construct_classic_plane_quadric_from_edge(typename boost::graph_traits<TriangleMesh>::halfedge_descriptor he,
const TriangleMesh& mesh,
const VertexPointMap vpm,
const GeomTraits& gt)
{
typedef typename GeomTraits::Vector_3 Vector_3;
typedef typename boost::graph_traits<TriangleMesh>::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 TriangleMesh, typename VertexPointMap, typename GeomTraits>
typename GarlandHeckbert_matrix_types<GeomTraits>::Mat_4
construct_classic_plane_quadric_from_face(typename boost::graph_traits<TriangleMesh>::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 TriangleMesh, typename VertexPointMap, typename GeomTraits>
typename GarlandHeckbert_matrix_types<GeomTraits>::Mat_4
construct_classic_triangle_quadric_from_face(typename boost::graph_traits<TriangleMesh>::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<GeomTraits>::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 GeomTraits>
typename GarlandHeckbert_matrix_types<GeomTraits>::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<GeomTraits>::Col_3 Col_3;
typedef typename GarlandHeckbert_matrix_types<GeomTraits>::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 TriangleMesh, typename VertexPointMap, typename GeomTraits>
typename GarlandHeckbert_matrix_types<GeomTraits>::Mat_4
construct_prob_triangle_quadric_from_face(typename boost::graph_traits<TriangleMesh>::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<GeomTraits>::Mat_3 Mat_3;
typedef typename GarlandHeckbert_matrix_types<GeomTraits>::Col_3 Col_3;
typedef typename GarlandHeckbert_matrix_types<GeomTraits>::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<GeomTraits>(a_minus_b);
const Mat_3 cp_bc = skew_sym_mat_cross_product<GeomTraits>(b_minus_c);
const Mat_3 cp_ca = skew_sym_mat_cross_product<GeomTraits>(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<GeomTraits>(cross_product(a_minus_b, ab))
+ vector_to_col_3<GeomTraits>(cross_product(b_minus_c, bc))
+ vector_to_col_3<GeomTraits>(cross_product(c_minus_a, ca)))
+ 2 * square(var) * vector_to_col_3<GeomTraits>(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 <typename TriangleMesh, typename GeomTraits>
std::pair<typename GeomTraits::FT, typename GeomTraits::FT>
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;

283
Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/internal/GarlandHeckbert_policy_base.h
View File

@ -28,149 +28,138 @@ namespace CGAL {
namespace Surface_mesh_simplification {
namespace internal {
template <typename Mat_4>
Mat_4 combine_matrices(const Mat_4& a, const Mat_4& b)
template <typename Mat>
Mat combine_matrices(const Mat& a, const Mat& b)
{
return a + b;
}
template <typename VertexCostMap,
typename GeomTraits,
typename QuadricImpl>
class GarlandHeckbert_cost_base
template <typename GeomTraits>
struct GarlandHeckbert_matrix_types
{
typedef typename GeomTraits::FT FT;
typedef Eigen::Matrix<FT, 3, 3, Eigen::DontAlign> Mat_3;
typedef Eigen::Matrix<FT, 3, 1> Col_3;
typedef Eigen::Matrix<FT, 4, 4, Eigen::DontAlign> Mat_4;
typedef Eigen::Matrix<FT, 4, 1> Col_4;
typedef Eigen::Matrix<FT, 1, 4> Row_4;
};
// Storage is initialized by the most-derived class (e.g. GarlandHeckbert_plane_policies)
template <typename QuadricCalculator, typename TriangleMesh, typename GeomTraits>
struct GarlandHeckbert_quadrics_storage
{
typedef typename GarlandHeckbert_matrix_types<GeomTraits>::Mat_4 Mat_4;
typedef typename GarlandHeckbert_matrix_types<GeomTraits>::Col_4 Col_4;
typedef Mat_4 Cost_matrix;
typedef CGAL::dynamic_vertex_property_t<Cost_matrix> Cost_property;
typedef typename boost::property_map<TriangleMesh, Cost_property>::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 <typename QuadricCalculator, typename TriangleMesh, typename GeomTraits>
class GarlandHeckbert_cost_base
// Virtual base because the storage is common to both the cost and the placement
: virtual public GarlandHeckbert_quadrics_storage<QuadricCalculator, TriangleMesh, GeomTraits>
{
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<TriangleMesh>::vertex_descriptor vertex_descriptor;
typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::graph_traits<TriangleMesh>::face_descriptor face_descriptor;
typedef typename GeomTraits::FT FT;
typedef Eigen::Matrix<FT, 4, 4, Eigen::DontAlign> Mat_4;
typedef Eigen::Matrix<FT, 4, 1> 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 <typename TM, typename VPM>
Mat_4 construct_quadric(typename boost::graph_traits<TM>::face_descriptor f,
const TM& tmesh,
const VPM point_map,
const GeomTraits& gt) const
{
return static_cast<const QuadricImpl*>(this)->construct_quadric_from_face(f, tmesh, point_map, gt);
}
template <typename TM, typename VPM>
Mat_4 construct_quadric(typename boost::graph_traits<TM>::halfedge_descriptor he,
const TM& tmesh,
const VPM point_map,
const GeomTraits& gt) const
{
return static_cast<const QuadricImpl*>(this)->construct_quadric_from_edge(he, tmesh, point_map, gt);
}
template <typename TM, typename VPM>
Vector_3 construct_edge_normal(typename boost::graph_traits<TM>::halfedge_descriptor he,
const TM& tmesh,
const VPM point_map,
const GeomTraits& gt) const
{
typedef typename boost::graph_traits<TM>::vertex_descriptor vertex_descriptor;
const Vector_3 face_normal = construct_unit_normal_from_face(face(he, tmesh), tmesh, point_map, gt);
const vertex_descriptor vs = source(he, tmesh);
const vertex_descriptor vt = target(he, tmesh);
const Vector_3 edge_vector { get(point_map, vs), get(point_map, vt) };
const Vector_3 discontinuity_normal = cross_product(edge_vector, face_normal);
// normalize
const Vector_3 normal = discontinuity_normal / sqrt(discontinuity_normal.squared_length());
return normal;
}
template <typename VPM, typename TM>
Vector_3 construct_unit_normal_from_face(typename boost::graph_traits<TM>::face_descriptor f,
const TM& tmesh,
const VPM point_map,
const GeomTraits& gt) const
{
typedef typename boost::graph_traits<TM>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::property_traits<VPM>::reference Point_reference;
auto unit_normal = gt.construct_unit_normal_3_object();
const halfedge_descriptor h = halfedge(f, tmesh);
// get the three points of the face and calculate their unit normal
const Point_reference p = get(point_map, source(h, tmesh));
const Point_reference q = get(point_map, target(h, tmesh));
const Point_reference r = get(point_map, target(next(h, tmesh), tmesh));
return unit_normal(p, q, r);
}
template <typename TM>
static bool is_discontinuity_edge(const typename boost::graph_traits<TM>::halfedge_descriptor h,
const TM& tmesh)
static bool is_discontinuity_edge(const halfedge_descriptor h,
const TriangleMesh& tmesh)
{
return is_border_edge(h, tmesh);
}
public:
template <typename VertexPointMap>
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 <typename VertexPointMap>
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 <typename TM, typename VPM>
void initialize(const TM& tmesh,
const VPM vpm,
template <typename VertexPointMap>
void initialize(const TriangleMesh& tmesh,
const VertexPointMap vpm,
const GeomTraits& gt) const
{
typedef boost::graph_traits<TM> 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<VPM>::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<TriangleMesh>::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 <typename Profile, typename VertexDescriptor>
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 <typename Profile>
boost::optional<typename Profile::FT>
operator()(const Profile& profile,
@ -205,88 +203,69 @@ public:
typedef boost::optional<typename Profile::FT> Optional_FT;
if(!placement)
{
// return empty value
return boost::optional<typename Profile::FT>();
}
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 <typename Profile, typename VertexDescriptor>
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 <typename VertexCostMap,
typename GeomTraits,
typename QuadricImpl>
template <typename QuadricCalculator, typename TriangleMesh, typename GeomTraits>
class GarlandHeckbert_placement_base
// Virtual base because the storage is common to both the cost and the placement
: virtual public GarlandHeckbert_quadrics_storage<QuadricCalculator, TriangleMesh, GeomTraits>
{
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<FT, 4, 4, Eigen::DontAlign> Mat_4;
typedef Eigen::Matrix<FT, 4, 1> Col_4;
private:
Vertex_cost_map m_cost_matrices;
typedef QuadricCalculator Quadric_calculator;
typedef typename GarlandHeckbert_matrix_types<GeomTraits>::Mat_4 Mat_4;
typedef typename GarlandHeckbert_matrix_types<GeomTraits>::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<const QuadricImpl*>(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 <typename Profile>
boost::optional<typename Profile::Point> operator()(const Profile& profile) const
{
CGAL_precondition(!get(m_cost_matrices, profile.v0()).isZero(0));
CGAL_precondition(!get(m_cost_matrices, profile.v1()).isZero(0));
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<typename Profile::Point> pt = typename Profile::Point(opt(0) / opt(3),
opt(1) / opt(3),