Make it work with Heat_method_traits_3

This commit is contained in:
Andreas Fabri 2018-10-30 22:31:24 +01:00
parent 347fab4d34
commit 470ce56df5
10 changed files with 280 additions and 52 deletions

View File

@ -32,7 +32,11 @@ public:
typedef unspecified_type Compute_y_3;
typedef unspecified_type Compute_z_3;
typedef unspecified_type Construct_vector_3;
typedef unspecified_type Construct_sum_of_vectors_3;
typedef unspecified_type Construct_scaled_vector_3;
typedef unspecified_type Construct_cross_product_vector_3;
typedef unspecified_type Compute_scalar_product_3;
typedef unspecified_type Compute_squared_distance_3;
/// @}

View File

@ -33,7 +33,8 @@ step will roughly double the overall preprocessing cost.
In the next section we give some examples. Section \ref sec_HM_definitions presents
the mathematical theory of the Heat method. The last section is about the \ref sec_HM_history.
Note that this package requires the third party \ref thirdpartyEigen library.
Note that this package depends on the third party \ref thirdpartyEigen library, or another
model of the concept `SparseLinearAlgebraWithFactorTraits_d`.
This implementation is based on \cgalCite{cgal:cww-ghnac-13} and \cgalCite{cgal:fsbs-acidt-06}

View File

@ -37,6 +37,8 @@ source vertices. }
## Classes ##
- `CGAL::Heat_method_3::Surface_mesh_geodesic_distances_3`
- `CGAL::Heat_method_3::Direct`
- `CGAL::Heat_method_3::Intrinsic_Delaunay`
## Functions ##

View File

@ -54,6 +54,7 @@ include_directories( BEFORE include )
include( CGAL_CreateSingleSourceCGALProgram )
create_single_source_cgal_program( "heat_method.cpp" )
create_single_source_cgal_program( "heat_method_concept.cpp" )
create_single_source_cgal_program( "heat_method_polyhedron.cpp" )
create_single_source_cgal_program( "heat_method_surface_mesh.cpp" )
create_single_source_cgal_program( "heat_method_surface_mesh_intrinsic.cpp" )

View File

@ -231,6 +231,7 @@ private:
double
summation_of_edges() const
{
typename Traits::Compute_squared_distance_3 squared_distance = Traits().compute_squared_distance_3_object();
double edge_sum = 0;
CGAL::Vertex_around_face_iterator<TriangleMesh> vbegin, vend, vmiddle;
BOOST_FOREACH(face_descriptor f, faces(tm)) {
@ -245,9 +246,9 @@ private:
VertexPointMap_reference pi = get(vpm,current);
VertexPointMap_reference pj = get(vpm, neighbor_one);
VertexPointMap_reference pk = get(vpm, neighbor_two);
edge_sum += (CGAL::is_border(opposite(hd,tm),tm)?1.0:0.5) * std::sqrt(CGAL::squared_distance(pi,pj));
edge_sum += (CGAL::is_border(opposite(hd2,tm),tm)?1.0:0.5) * std::sqrt(CGAL::squared_distance(pj,pk)) ;
edge_sum += (CGAL::is_border(opposite(hd3,tm),tm)?1.0:0.5) * std::sqrt(CGAL::squared_distance(pk,pi)) ;
edge_sum += (CGAL::is_border(opposite(hd,tm),tm)?1.0:0.5) * std::sqrt(squared_distance(pi,pj));
edge_sum += (CGAL::is_border(opposite(hd2,tm),tm)?1.0:0.5) * std::sqrt(squared_distance(pj,pk)) ;
edge_sum += (CGAL::is_border(opposite(hd3,tm),tm)?1.0:0.5) * std::sqrt(squared_distance(pk,pi)) ;
}
return edge_sum;
}
@ -312,6 +313,11 @@ private:
void
compute_unit_gradient()
{
Traits::Construct_vector_3 construct_vector = Traits().construct_vector_3_object();
Traits::Construct_sum_of_vectors_3 sum = Traits().construct_sum_of_vectors_3_object();
Traits::Compute_scalar_product_3 scalar_product = Traits().compute_scalar_product_3_object();
Traits::Construct_cross_product_vector_3 cross_product = Traits().construct_cross_product_vector_3_object();
Traits::Construct_scaled_vector_3 scale = Traits().construct_scaled_vector_3_object();
if(X.empty()){
X.resize(num_faces(tm));
}
@ -335,9 +341,13 @@ private:
//cross that with eij, ejk, eki
//so (Ncross eij) *uk and so on
//sum all of those then multiply by 1./(2a)
Vector_3 cross = CGAL::cross_product((p_j-p_i), (p_k-p_i));
double N_cross = (CGAL::sqrt(cross*cross));
Vector_3 unit_cross = cross/N_cross;
Vector_3 v_ij = construct_vector(p_i,p_j);
Vector_3 v_ik = construct_vector(p_i,p_k);
Vector_3 cross = cross_product(v_ij, v_ik);
double N_cross = (CGAL::sqrt(scalar_product(cross,cross)));
Vector_3 unit_cross = scale(cross, 1./N_cross);
double area_face = N_cross * (1./2);
double u_i = CGAL::abs(solved_u(i));
double u_j = CGAL::abs(solved_u(j));
@ -349,25 +359,24 @@ private:
u_j = u_j * r_Mag;
u_k = u_k * r_Mag;
}
Vector_3 edge_sums = u_k * CGAL::cross_product(unit_cross,(p_j-p_i));
edge_sums = edge_sums + u_i * (CGAL::cross_product(unit_cross, (p_k-p_j)));
edge_sums = edge_sums + u_j * CGAL::cross_product(unit_cross, (p_i-p_k));
edge_sums = edge_sums * (1./area_face);
double e_magnitude = CGAL::sqrt(edge_sums*edge_sums);
X[face_i] = edge_sums*(1./e_magnitude);
Vector_3 edge_sums = scale(cross_product(unit_cross,v_ij), u_k);
edge_sums = sum(edge_sums, scale(cross_product(unit_cross, construct_vector(p_j,p_k)), u_i));
edge_sums = sum(edge_sums, scale(cross_product(unit_cross, construct_vector(p_k,p_i)), u_j));
edge_sums = scale(edge_sums, (1./area_face));
double e_magnitude = CGAL::sqrt(scalar_product(edge_sums,edge_sums));
X[face_i] = scale(edge_sums,(1./e_magnitude));
}
}
double
dot_eigen_vector(const Vector_3& a, const Vector_3& b) const
{
return (a.x()*b.x() + a.y()*b.y() + a.z()*b.z());
}
void
compute_divergence()
{
Traits::Construct_cross_product_vector_3 cross_product = Traits().construct_cross_product_vector_3_object();
Traits::Compute_scalar_product_3 scalar_product = Traits().compute_scalar_product_3_object();
Traits::Construct_vector_3 construct_vector = Traits().construct_vector_3_object();
Matrix indexD(dimension,1);
CGAL::Vertex_around_face_iterator<TriangleMesh> vbegin, vend, vmiddle;
BOOST_FOREACH(face_descriptor f, faces(tm)) {
@ -382,24 +391,32 @@ private:
VertexPointMap_reference p_j = get(vpm, neighbor_one);
VertexPointMap_reference p_k = get(vpm, neighbor_two);
Index face_i = get(face_id_map, f);
Vector_3 cross = CGAL::cross_product((p_j-p_i), (p_k-p_i));
double norm_cross = (CGAL::sqrt(cross*cross));
double dot = (p_j-p_i)*(p_k-p_i);
Vector_3 v_ij = construct_vector(p_i,p_j);
Vector_3 v_ik = construct_vector(p_i,p_k);
Vector_3 cross = cross_product(v_ij, v_ik);
double norm_cross = CGAL::sqrt(scalar_product(cross,cross));
double dot = scalar_product(v_ij, v_ik);
double cotan_i = dot/norm_cross;
cross = CGAL::cross_product((p_i-p_j), (p_k-p_j));
dot = to_double((p_i-p_j)*(p_k-p_j));
Vector_3 v_ji = construct_vector(p_j, p_i);
Vector_3 v_jk = construct_vector(p_j, p_k);
cross = cross_product(v_ji, v_jk);
dot = to_double(scalar_product(v_ji, v_jk));
double cotan_j = dot/norm_cross;
cross = CGAL::cross_product((p_i-p_k), (p_j-p_k));
dot = to_double((p_i-p_k)*(p_j-p_k));
Vector_3 v_ki = construct_vector(p_k,p_i);
Vector_3 v_kj = construct_vector(p_k,p_j);
cross = cross_product(v_ki, v_kj);
dot = to_double(scalar_product(v_ki,v_kj));
double cotan_k = dot/norm_cross;
const Vector_3& a = X[face_i];
double i_entry = cotan_k*(dot_eigen_vector(a,(p_j-p_i))) + cotan_j*(dot_eigen_vector(a,(p_k-p_i)));
double j_entry = cotan_i*(dot_eigen_vector(a,(p_k-p_j))) + cotan_k*(dot_eigen_vector(a,(p_i-p_j)));
double k_entry = cotan_j*(dot_eigen_vector(a,(p_i-p_k))) + cotan_i*(dot_eigen_vector(a,(p_j-p_k)));
double i_entry = (scalar_product(a,v_ij) * cotan_k) + (scalar_product(a,v_ik) * cotan_j);
double j_entry = (scalar_product(a,v_jk) * cotan_i) + (scalar_product(a,v_ji) * cotan_k);
double k_entry = (scalar_product(a,v_ki) * cotan_j) + (scalar_product(a,v_kj) * cotan_i);
indexD.add_coef(i, 0, (1./2)*i_entry);
indexD.add_coef(j, 0, (1./2)*j_entry);
@ -499,6 +516,10 @@ private:
void
build()
{
Traits::Construct_cross_product_vector_3 cross_product = Traits().construct_cross_product_vector_3_object();
Traits::Compute_scalar_product_3 scalar_product = Traits().compute_scalar_product_3_object();
Traits::Construct_vector_3 construct_vector = Traits().construct_vector_3_object();
source_change_flag = false;
CGAL_precondition(is_triangle_mesh(tm));
@ -540,27 +561,35 @@ private:
pj = p_j;
pk = p_k;
Vector_3 cross = CGAL::cross_product((pj-pi), (pk-pi));
double dot = (pj-pi)*(pk-pi);
double norm_cross = (CGAL::sqrt(cross*cross));
Vector_3 v_ij = construct_vector(p_i,p_j);
Vector_3 v_ik = construct_vector(p_i,p_k);
Vector_3 cross = cross_product(v_ij, v_ik);
double dot = scalar_product(v_ij,v_ik);
double norm_cross = (CGAL::sqrt(scalar_product(cross,cross)));
double cotan_i = dot/norm_cross;
m_cotan_matrix.add_coef(j,k ,-(1./2)*cotan_i);
m_cotan_matrix.add_coef(k,j,-(1./2)* cotan_i);
m_cotan_matrix.add_coef(j,j,(1./2)*cotan_i);
m_cotan_matrix.add_coef(k,k,(1./2)* cotan_i);
cross = CGAL::cross_product((pi-pj), (pk-pj));
dot = to_double((pi-pj)*(pk-pj));
Vector_3 v_ji = construct_vector(p_j,p_i);
Vector_3 v_jk = construct_vector(p_j,p_k);
cross = cross_product(v_ji, v_jk);
dot = to_double(scalar_product(v_ji, v_jk));
double cotan_j = dot/norm_cross;
m_cotan_matrix.add_coef(i,k ,-(1./2)*cotan_j);
m_cotan_matrix.add_coef(k,i,-(1./2)* cotan_j);
m_cotan_matrix.add_coef(i,i,(1./2)* cotan_j);
m_cotan_matrix.add_coef(k,k,(1./2)* cotan_j);
cross = CGAL::cross_product((pi-pk), (pj-pk));
dot = to_double((pi-pk)*(pj-pk));
Vector_3 v_ki = construct_vector(p_k,p_i);
Vector_3 v_kj = construct_vector(p_k,p_j);
cross = cross_product(v_ki, v_kj);
dot = to_double(scalar_product(v_ki,v_kj));
double cotan_k = dot/norm_cross;
m_cotan_matrix.add_coef(i,j,-(1./2)*cotan_k);
m_cotan_matrix.add_coef(j,i,-(1./2)* cotan_k);
@ -604,9 +633,7 @@ private:
Vector solved_phi;
std::set<Index> source_index;
bool source_change_flag;
std::vector<Eigen::Triplet<double> > m_triplets;
};
template <typename TriangleMesh,
@ -859,7 +886,8 @@ public:
/// \ingroup PkgHeatMethod
/// computes for each vertex of the triangle mesh `tm` the estimated geodesic distance to a given source vertex.
/// \tparam TriangleMesh a triangulated surface mesh, model of `FaceListGraph` and `HalfedgeListGraph`
/// \tparam TriangleMesh a triangulated surface mesh, model of `FaceListGraph` and `HalfedgeListGraph`.
/// It must have an internal vertex point property map with the value type being a 3D point from a cgal Kernel model
/// \tparam VertexDistanceMap a property map model of `WritablePropertyMap`
/// with `boost::graph_traits<TriangleMesh>::vertex_descriptor` as key type and `double` as value type.
/// \tparam Mode either the tag `Direct` or `Intrinsic_Delaunay`, which determines if the geodesic distance
@ -867,7 +895,7 @@ public:
/// The default is `Direct`.
///
/// \sa CGAL::Heat_method_3::Surface_mesh_geodesic_distances_3
template <typename TriangleMesh, typename VertexDistanceMap,Mode>
template <typename TriangleMesh, typename VertexDistanceMap, typename Mode>
void
estimate_geodesic_distances(const TriangleMesh& tm,
VertexDistanceMap vdm,
@ -879,7 +907,8 @@ estimate_geodesic_distances(const TriangleMesh& tm,
hm.estimate_geodesic_distances(vdm);
}
#ifndef DOXYGEN_RUNNING
template <typename TriangleMesh, typename VertexDistanceMap>
void
estimate_geodesic_distances(const TriangleMesh& tm,
@ -890,11 +919,13 @@ estimate_geodesic_distances(const TriangleMesh& tm,
hm.add_source(source);
hm.estimate_geodesic_distances(vdm);
}
#endif
/// \ingroup PkgHeatMethod
/// \ingroup PkgHeatMethod
/// computes for each vertex of the triangle mesh `tm` the estimated geodesic distance to a given set of source vertices.
/// \tparam TriangleMesh a triangulated surface mesh, model of `FaceListGraph` and `HalfedgeListGraph`
/// It must have an internal vertex point property map with the value type being a 3D point from a cgal Kernel model
/// \tparam VertexDistanceMap a property map model of `WritablePropertyMap`
/// with `boost::graph_traits<TriangleMesh>::vertex_descriptor` as key type and `double` as value type.
/// \tparam VertexRange a range with value type `boost::graph_traits<TriangleMesh>::vertex_descriptor`

View File

@ -24,6 +24,8 @@
#include <CGAL/license/Heat_method_3.h>
#ifndef DOXYGEN_RUNNING
namespace CGAL {
namespace Heat_method_3 {
@ -39,7 +41,7 @@ namespace Heat_method_3 {
return t;
}
};
#endif
} // namespace Heat_method_3
} // namespace CGAL
#endif //CGAL_HEAT_METHOD_3_INTERNAL_V2V_H

View File

@ -53,6 +53,7 @@ include( CGAL_CreateSingleSourceCGALProgram )
create_single_source_cgal_program( "validate.cpp" )
create_single_source_cgal_program( "idt.cpp" )
create_single_source_cgal_program( "heat_method_concept.cpp" )
create_single_source_cgal_program( "heat_method_surface_mesh_test.cpp" )
create_single_source_cgal_program( "heat_method_intrinsic_surface_mesh_test.cpp" )
create_single_source_cgal_program( "intrinsic_delaunay_triangulation_test.cpp" )

View File

@ -0,0 +1,123 @@
struct Heat_method_traits_3
{
typedef double FT;
struct Point_3 {
Point_3()
{}
Point_3(double, double,double)
{}
Point_3(const Point_3&)
{}
};
struct Vector_3 {
Vector_3()
{}
Vector_3(double, double, double)
{}
Vector_3(const Vector_3&)
{}
};
struct Point_2 {
Point_2()
{}
Point_2(double, double)
{}
Point_2(const Point_2&)
{}
};
struct Compute_x_2 {
double operator()(const Point_2&)const
{return 0;}
};
struct Compute_y_2 {
double operator()(const Point_2&)const
{return 0;}
};
struct Construct_vector_3{
Vector_3 operator()(const Point_3&, const Point_3&) const
{ return Vector_3();}
};
struct Construct_scaled_vector_3{
Vector_3 operator()(const Vector_3&, double) const
{ return Vector_3();}
};
struct Construct_cross_product_vector_3 {
Vector_3 operator()(const Vector_3&, const Vector_3&) const
{ return Vector_3();}
};
struct Construct_sum_of_vectors_3 {
Vector_3 operator()(const Vector_3&, const Vector_3&) const
{ return Vector_3();}
};
struct Compute_scalar_product_3 {
double operator()(const Vector_3&, const Vector_3&) const
{ return 0;}
};
struct Compute_squared_distance_3 {
double operator()(const Point_3&, const Point_3&) const
{ return 0;}
};
Compute_x_2 compute_x_2_object() const
{
return Compute_x_2();
}
Compute_y_2 compute_y_2_object() const
{
return Compute_y_2();
}
Construct_vector_3 construct_vector_3_object() const
{
return Construct_vector_3();
}
Construct_scaled_vector_3 construct_scaled_vector_3_object() const
{
return Construct_scaled_vector_3();
}
Construct_cross_product_vector_3 construct_cross_product_vector_3_object() const
{
return Construct_cross_product_vector_3();
}
Compute_squared_distance_3 compute_squared_distance_3_object() const
{
return Compute_squared_distance_3();
}
Compute_scalar_product_3 compute_scalar_product_3_object() const
{
return Compute_scalar_product_3();
}
Construct_sum_of_vectors_3 construct_sum_of_vectors_3_object() const
{
return Construct_sum_of_vectors_3();
}
};

View File

@ -0,0 +1,63 @@
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Heat_method_3/Surface_mesh_geodesic_distances_3.h>
#include "Heat_method_3/Heat_method_traits_3.h"
#include <fstream>
#include <iostream>
#include <boost/foreach.hpp>
typedef Heat_method_traits_3 Kernel;
typedef Kernel::Point_3 Point_3;
typedef CGAL::Surface_mesh<Point_3> Surface_mesh;
typedef boost::graph_traits<Surface_mesh>::vertex_descriptor vertex_descriptor;
typedef Surface_mesh::Property_map<vertex_descriptor,double> Vertex_distance_map;
typedef CGAL::Heat_method_3::Surface_mesh_geodesic_distances_3<Surface_mesh,
CGAL::Heat_method_3::Direct,
CGAL::Eigen_solver_traits<Eigen::SimplicialLDLT<CGAL::Eigen_sparse_matrix<double>::EigenType > >,
boost::property_map< Surface_mesh, CGAL::vertex_point_t>::const_type,
Kernel> Heat_method;
int main(int argc, char* argv[])
{
//read in mesh
Surface_mesh sm;
//the heat intensity will hold the distance values from the source set
Vertex_distance_map heat_intensity = sm.add_property_map<vertex_descriptor, double>("v:heat_intensity", 0).first;
Heat_method hm(sm);
//add the first vertex as the source set
vertex_descriptor source = *(vertices(sm).first);
hm.add_source(source);
hm.estimate_geodesic_distances(heat_intensity);
Point_3 sp = sm.point(source);
vertex_descriptor far;
double sdistance = 0;
BOOST_FOREACH(vertex_descriptor vd , vertices(sm)){
std::cout << vd << " is at distance " << get(heat_intensity, vd) << " from " << source << std::endl;
/*
if(squared_distance(sp,sm.point(vd)) > sdistance){
far = vd;
sdistance = squared_distance(sp,sm.point(vd));
}
*/
}
hm.add_source(far);
hm.estimate_geodesic_distances(heat_intensity);
BOOST_FOREACH(vertex_descriptor vd , vertices(sm)){
std::cout << vd << " is at distance " << get(heat_intensity, vd) << " " << std::endl;
}
return 0;
}

View File

@ -10,7 +10,7 @@
#include <QMessageBox>
#include <CGAL/boost/graph/helpers.h>
#include <CGAL/Polygon_mesh_processing/compute_normal.h>
#include <CGAL/Heat_method_3/Heat_method_3.h>
#include <CGAL/Heat_method_3/Surface_mesh_geodesic_distances_3.h>
#include "Scene_points_with_normal_item.h"
#include "Messages_interface.h"
@ -309,8 +309,8 @@ class DisplayPropertyPlugin :
Q_INTERFACES(CGAL::Three::Polyhedron_demo_plugin_interface)
Q_PLUGIN_METADATA(IID "com.geometryfactory.PolyhedronDemo.PluginInterface/1.0")
typedef SMesh::Property_map<boost::graph_traits<SMesh>::vertex_descriptor, double> Vertex_distance_map;
typedef CGAL::Heat_method_3::Heat_method_3<SMesh> Heat_method;
typedef CGAL::Heat_method_3::Heat_method_3<SMesh, CGAL::Tag_true> Heat_method_idt;
typedef CGAL::Heat_method_3::Surface_mesh_geodesic_distances_3<SMesh> Heat_method;
typedef CGAL::Heat_method_3::Surface_mesh_geodesic_distances_3<SMesh, CGAL::Heat_method_3::Intrinsic_Delaunay> Heat_method_idt;
public:
@ -769,9 +769,9 @@ private Q_SLOTS:
}
if(iDT){
hm_idt->fill_distance_map(heat_intensity);
hm_idt->estimate_geodesic_distances(heat_intensity);
}else{
hm->fill_distance_map(heat_intensity);
hm->estimate_geodesic_distances(heat_intensity);
}
double max = 0;