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

Improve reader sanity 

- Re-indent to more normal standards
- Remove usage of *_impl.h
- Clean useless and trailing whitespace
- some minor code readability changes
This commit is contained in:
Mael Rouxel-Labbé 2019-05-03 14:56:21 +02:00
parent c8a3fdaff6
commit 4747f8ccaf
54 changed files with 3052 additions and 3443 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Surface_mesh_simplification/doc/Surface_mesh_simplification/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Bounded_normal_change_placement.h
View File

@ -44,7 +44,7 @@ triangle in the profile changes the normal by more than 90 degree.
*/
template <typename Profile>
optional<typename Profile::Point>
operator()( Profile const& profile ) const;
operator()(const Profile& profile) const;
/// @}

2
Surface_mesh_simplification/doc/Surface_mesh_simplification/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Constrained_placement.h
View File

@ -30,7 +30,7 @@ Constructor
*/
Constrained_placement(
EdgeIsConstrainedMap map=EdgeIsConstrainedMap(),
BasePlacement base= BasePlacement() );
BasePlacement base= BasePlacement());
/// @}
}; /* end Surface_mesh_simplification::Midpoint_placement */

10
Surface_mesh_simplification/doc/Surface_mesh_simplification/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Count_ratio_stop_predicate.h
View File

@ -26,7 +26,7 @@ public:
/*!
Initializes the predicate establishing the `ratio`.
*/
Count_ratio_stop_predicate<TriangleMesh>( double ratio );
Count_ratio_stop_predicate<TriangleMesh>(double ratio);
/// @}
@ -34,14 +34,14 @@ Count_ratio_stop_predicate<TriangleMesh>( double ratio );
/// @{
/*!
Returns ` ( ((double)current_count / (double)initial_count) < ratio)`.
Returns ` (((double)current_count / (double)initial_count) < ratio)`.
All other parameters are ignored (but exist since this is a generic policy).
*/
bool operator()( FT const& current_cost
, Profile const& edge_profile
bool operator()(FT const& current_cost
, const Profile& edge_profile
, size_type initial_count
, size_type current_count
) const ;
) const;
/// @}

8
Surface_mesh_simplification/doc/Surface_mesh_simplification/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Count_stop_predicate.h
View File

@ -25,7 +25,7 @@ public:
/*!
Initializes the predicate establishing the `threshold` value.
*/
Count_stop_predicate<TriangleMesh>( size_type threshold );
Count_stop_predicate<TriangleMesh>(size_type threshold);
/// @}
@ -35,11 +35,11 @@ Count_stop_predicate<TriangleMesh>( size_type threshold );
/*!
Returns `(current_count < threshold)`. All other parameters are ignored (but exist since this is a generic policy).
*/
bool operator()( FT const& current_cost
, Profile const& edge_profile
bool operator()(FT const& current_cost
, const Profile& edge_profile
, size_type initial_count
, size_type current_count
) const ;
) const;
/// @}

4
Surface_mesh_simplification/doc/Surface_mesh_simplification/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_length_cost.h
View File

@ -39,8 +39,8 @@ The `placement` argument is ignored.
*/
template <typename Profile, typename T>
optional<typename Profile::FT> operator()( Profile const& profile
, T const& placement ) const;
optional<typename Profile::FT> operator()(const Profile& profile
, T const& placement) const;
/// @}

8
Surface_mesh_simplification/doc/Surface_mesh_simplification/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_length_stop_predicate.h
View File

@ -24,7 +24,7 @@ public:
/*!
Initializes the predicate establishing the `threshold` value.
*/
Edge_length_stop_predicate<TriangleMesh>( FT threshold );
Edge_length_stop_predicate<TriangleMesh>(FT threshold);
/// @}
@ -35,11 +35,11 @@ Edge_length_stop_predicate<TriangleMesh>( FT threshold );
Returns `(CGAL::squared_distance(edge_profile.p0(),edge_profile.p1()) > threshold*threshold)`.
All other parameters are ignored (but exist since this is a generic policy).
*/
bool operator()( FT const&
, Profile const& edge_profile
bool operator()(FT const&
, const Profile& edge_profile
, size_type
, size_type
) const ;
) const;
/// @}

6
Surface_mesh_simplification/doc/Surface_mesh_simplification/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/LindstromTurk_cost.h
View File

@ -26,7 +26,7 @@ public:
Initializes the policy with the given <I>weighting unit factor</I>.
See \ref SurfaceMeshSimplificationLindstromTurkStrategy for details on the meaning of this factor.
*/
LindstromTurk_cost<TriangleMesh>( FT const& factor = FT(0.5) );
LindstromTurk_cost<TriangleMesh>(FT const& factor = FT(0.5));
/// @}
@ -39,8 +39,8 @@ the new `placement` computed for it.
*/
template <typename Profile>
optional<typename Profile::FT>
operator()( Profile const& profile
, boost::optional<typename Profile::Point> const& placement ) const;
operator()(const Profile& profile
, boost::optional<typename Profile::Point> const& placement) const;
/// @}

4
Surface_mesh_simplification/doc/Surface_mesh_simplification/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/LindstromTurk_placement.h
View File

@ -28,7 +28,7 @@ public:
Initializes the policy with the given <I>weighting unit factor</I>.
See \ref SurfaceMeshSimplificationLindstromTurkStrategy for details on the meaning of this factor.
*/
LindstromTurk_placement<TriangleMesh>( FT const& factor = FT(0.5) );
LindstromTurk_placement<TriangleMesh>(FT const& factor = FT(0.5));
/// @}
@ -41,7 +41,7 @@ Returns the new position for the remaining vertex after collapsing the edge
*/
template <typename Profile>
optional<typename Profile::Point>
operator()( Profile const& profile ) const;
operator()(const Profile& profile) const;
/// @}

2
Surface_mesh_simplification/doc/Surface_mesh_simplification/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Midpoint_placement.h
View File

@ -37,7 +37,7 @@ the points of the source and target vertices
*/
template <typename Profile>
optional<typename Profile::Point>
operator()( Profile const& profile ) const;
operator()(const Profile& profile) const;
/// @}

16
Surface_mesh_simplification/doc/Surface_mesh_simplification/Concepts/EdgeCollapseSimplificationVisitor.h
View File

@ -49,26 +49,26 @@ typedef unspecified_type size_type;
/*!
Called before the algorithm starts.
*/
void OnStarted( TriangleMesh& surface_mesh );
void OnStarted(TriangleMesh& surface_mesh);
/*!
Called after the algorithm finishes.
*/
void OnFinished ( TriangleMesh& surface_mesh ) ;
void OnFinished (TriangleMesh& surface_mesh);
/*!
Called when the `StopPredicate` returned `true`
(but not if the algorithm terminates because the surface mesh could not be simplified any further).
*/
void OnStopConditionReached( TriangleMesh& surface_mesh ) ;
void OnStopConditionReached(TriangleMesh& surface_mesh);
/*!
Called during the <I>collecting phase</I> (when a cost is assigned to the edges),
for each edge collected.
*/
void OnCollected( Profile const& profile, boost::optional<FT> cost );
void OnCollected(const Profile& profile, boost::optional<FT> cost);
/*!
Called during the <I>processing phase</I> (when edges are collapsed),
@ -85,7 +85,7 @@ the edge will not be collapsed.
the number of edges.
*/
void OnSelected( Profile const& profile
void OnSelected(const Profile& profile
, boost::optional<FT> cost
, size_type initial_count
, size_type current_count
@ -99,14 +99,14 @@ If `placement` is absent (meaning that it could not be computed)
the edge will not be collapsed.
*/
void OnCollapsing( Profile const& profile
void OnCollapsing(const Profile& profile
, boost::optional<Point> placement
);
/*!
Called when an edge has been collapsed and replaced by the vertex `vd`
*/
void OnCollapsed( Profile const&, Profile::vertex_descriptor const& vd) {}
void OnCollapsed(Profile const&, Profile::const vertex_descriptor vd) {}
/*!
Called for each selected edge which cannot be
@ -115,7 +115,7 @@ type of the surface mesh (turn it into a non-manifold
for instance).
*/
void OnNonCollapsable( Profile const& profile );
void OnNonCollapsable(const Profile& profile);
/// @}

4
Surface_mesh_simplification/doc/Surface_mesh_simplification/Concepts/GetCost.h
View File

@ -32,8 +32,8 @@ using the calculated placement.
\tparam Profile must be a model of `EdgeProfile`.
*/
template <class Profile>
boost::optional<typename Profile::FT> operator()( Profile const& edge_profile
, boost::optional<typename Profile::Point> const& placement ) const;
boost::optional<typename Profile::FT> operator()(const Profile& edge_profile
, boost::optional<typename Profile::Point> const& placement) const;
/// @}

2
Surface_mesh_simplification/doc/Surface_mesh_simplification/Concepts/GetPlacement.h
View File

@ -34,7 +34,7 @@ which replaces the collapsing edge (represented by its profile).
\tparam Profile must be a model of `EdgeProfile`.
*/
template <class Profile>
boost::optional<Point>operator()( Profile const& edge_profile ) const;
boost::optional<Point>operator()(const Profile& edge_profile) const;
/// @}

6
Surface_mesh_simplification/doc/Surface_mesh_simplification/Concepts/StopPredicate.h
View File

@ -55,11 +55,11 @@ If the return value is `true` the simplification terminates before processing th
otherwise it continues normally.
*/
bool operator()( FT const& current_cost
, Profile const& profile
bool operator()(FT const& current_cost
, const Profile& profile
, size_type initial_count
, size_type current_count
) const ;
) const;
/// @}

10
Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_OpenMesh.cpp
View File

@ -53,18 +53,18 @@ private:
namespace SMS = CGAL::Surface_mesh_simplification ;
namespace SMS = CGAL::Surface_mesh_simplification;
int main( int argc, char** argv )
int main(int argc, char** argv)
{
Surface_mesh surface_mesh;
Constrained_edge_map constraints_map(surface_mesh);
if (argc==2)
if(argc==2)
OpenMesh::IO::read_mesh(surface_mesh, argv[1]);
else
OpenMesh::IO::read_mesh(surface_mesh, "cube.off");
if (!CGAL::is_triangle_mesh(surface_mesh)){
if(!CGAL::is_triangle_mesh(surface_mesh)){
std::cerr << "Input geometry is not triangulated." << std::endl;
return EXIT_FAILURE;
}
@ -93,7 +93,7 @@ int main( int argc, char** argv )
surface_mesh.garbage_collection();
std::cout << "\nFinished...\n" << r << " edges removed.\n"
<< num_edges(surface_mesh) << " final edges.\n" ;
<< num_edges(surface_mesh) << " final edges.\n";
OpenMesh::IO::write_mesh(surface_mesh, "out.off");

33
Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_all_short_edges.cpp
View File

@ -12,15 +12,14 @@
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_length_stop_predicate.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Midpoint_placement.h>
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef CGAL::Polyhedron_3<Kernel> Surface_mesh;
namespace SMS = CGAL::Surface_mesh_simplification ;
namespace SMS = CGAL::Surface_mesh_simplification;
int main( int argc, char** argv )
int main(int argc, char** argv)
{
if (argc<3)
if(argc<3)
{
std::cerr << "Usage: " << argv[0] << " input.off minimal_edge_length [out.off]\n";
return EXIT_FAILURE;
@ -28,29 +27,27 @@ int main( int argc, char** argv )
Surface_mesh surface_mesh;
std::ifstream is(argv[1]) ; is >> surface_mesh ;
std::ifstream is(argv[1]); is >> surface_mesh;
double threshold = atof(argv[2]);
if (!CGAL::is_triangle_mesh(surface_mesh)){
if(!CGAL::is_triangle_mesh(surface_mesh)){
std::cerr << "Input geometry is not triangulated." << std::endl;
return EXIT_FAILURE;
}
int r = SMS::edge_collapse
(surface_mesh
, CGAL::Surface_mesh_simplification::Edge_length_stop_predicate<double>(threshold)
, CGAL::parameters::vertex_index_map(get(CGAL::vertex_external_index,surface_mesh))
.halfedge_index_map (get(CGAL::halfedge_external_index ,surface_mesh))
.get_cost (SMS::Edge_length_cost <Surface_mesh>())
.get_placement(SMS::Midpoint_placement<Surface_mesh>())
);
int r = SMS::edge_collapse(surface_mesh,
CGAL::Surface_mesh_simplification::Edge_length_stop_predicate<double>(threshold),
CGAL::parameters::vertex_index_map(get(CGAL::vertex_external_index, surface_mesh))
.halfedge_index_map(get(CGAL::halfedge_external_index ,surface_mesh))
.get_cost(SMS::Edge_length_cost <Surface_mesh>())
.get_placement(SMS::Midpoint_placement<Surface_mesh>()));
std::cout << "\nFinished...\n" << r << " edges removed.\n"
<< (surface_mesh.size_of_halfedges()/2) << " final edges.\n" ;
<< (surface_mesh.size_of_halfedges()/2) << " final edges.\n";
std::ofstream os( argc > 3 ? argv[3] : "out.off" ) ;
std::ofstream os(argc > 3 ? argv[3] : "out.off");
os.precision(17);
os << surface_mesh ;
os << surface_mesh;
return EXIT_SUCCESS ;
return EXIT_SUCCESS;
}

12
Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_bounded_normal_change.cpp
View File

@ -9,13 +9,12 @@
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/LindstromTurk_placement.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Bounded_normal_change_placement.h>
typedef CGAL::Simple_cartesian<double> Kernel;
typedef CGAL::Surface_mesh<Kernel::Point_3> Surface_mesh;
namespace SMS = CGAL::Surface_mesh_simplification;
int main( int argc, char** argv )
int main(int argc, char** argv)
{
Surface_mesh surface_mesh;
@ -34,13 +33,12 @@ int main( int argc, char** argv )
// The surface mesh and stop conditions are mandatory arguments.
// The index maps are needed because the vertices and edges
// of this surface mesh lack an "id()" field.
SMS::edge_collapse( surface_mesh,
SMS::edge_collapse(surface_mesh,
stop,
CGAL::parameters::get_cost (SMS::LindstromTurk_cost<Surface_mesh>())
.get_placement(Placement())
);
CGAL::parameters::get_cost(SMS::LindstromTurk_cost<Surface_mesh>())
.get_placement(Placement()));
std::ofstream os( argc > 2 ? argv[2] : "out.off" );
std::ofstream os(argc > 2 ? argv[2] : "out.off");
os.precision(17);
os << surface_mesh;

123
Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_constrain_sharp_edges.cpp
View File

@ -1,6 +1,3 @@
#include <iostream>
#include <fstream>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Polyhedron_3.h>
@ -10,21 +7,25 @@
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Count_stop_predicate.h>
#include <CGAL/Unique_hash_map.h>
#include <CGAL/property_map.h>
#include <cmath>
#include <iostream>
#include <fstream>
typedef CGAL::Simple_cartesian<double> Kernel;
typedef Kernel::Point_3 Point_3;
typedef CGAL::Polyhedron_3<Kernel> Surface_mesh;
typedef boost::graph_traits<Surface_mesh>::vertex_descriptor vertex_descriptor;
typedef boost::graph_traits<Surface_mesh>::halfedge_descriptor halfedge_descriptor;
typedef boost::graph_traits<Surface_mesh>::edge_descriptor edge_descriptor;
typedef boost::graph_traits<Surface_mesh>::edge_iterator edge_iterator;
namespace SMS = CGAL::Surface_mesh_simplification ;
namespace SMS = CGAL::Surface_mesh_simplification;
//
// BGL property map which indicates whether an edge is marked as non-removable
struct Constrained_edge_map : public boost::put_get_helper<bool,Constrained_edge_map>
struct Constrained_edge_map
: public boost::put_get_helper<bool, Constrained_edge_map>
{
typedef boost::readable_property_map_tag category;
typedef bool value_type;
@ -37,7 +38,7 @@ struct Constrained_edge_map : public boost::put_get_helper<bool,Constrained_edge
reference operator[](key_type const& e) const { return is_constrained(e); }
bool is_constrained( key_type const& e ) const {
bool is_constrained(key_type const& e) const {
return mConstraints.is_defined(e);
}
@ -47,8 +48,8 @@ private:
bool is_border (edge_descriptor e, const Surface_mesh& sm)
{
return (face(halfedge(e,sm),sm) == boost::graph_traits<Surface_mesh>::null_face() )
|| (face(opposite(halfedge(e,sm),sm),sm) == boost::graph_traits<Surface_mesh>::null_face() );
return (face(halfedge(e,sm),sm) == boost::graph_traits<Surface_mesh>::null_face())
|| (face(opposite(halfedge(e,sm),sm),sm) == boost::graph_traits<Surface_mesh>::null_face());
}
Point_3 point(vertex_descriptor vd, const Surface_mesh& sm)
@ -56,26 +57,29 @@ Point_3 point(vertex_descriptor vd, const Surface_mesh& sm)
return get(CGAL::vertex_point, sm, vd);
}
int main( int argc, char** argv )
int main(int argc, char** argv)
{
CGAL::Unique_hash_map<edge_descriptor,bool> constraint_hmap(false);
Surface_mesh surface_mesh;
if (argc < 2){
if(argc < 2)
{
std::cerr << "Usage: " << argv[0] << " input.off [out.off]\n";
return EXIT_FAILURE;
}
std::ifstream is(argv[1]);
if(!is){
if(!is)
{
std::cerr << "Filename provided is invalid\n";
return EXIT_FAILURE;
}
is >> surface_mesh ;
if (!CGAL::is_triangle_mesh(surface_mesh)){
if(!CGAL::is_triangle_mesh(surface_mesh))
{
std::cerr << "Input geometry is not triangulated." << std::endl;
return EXIT_FAILURE;
}
@ -94,28 +98,29 @@ int main( int argc, char** argv )
// detect sharp edges
std::ofstream cst_output("constrained_edges.cgal");
edge_iterator eb,ee;
for(boost::tie(eb,ee) = edges(surface_mesh); eb != ee ; ++eb )
for(edge_descriptor ed : edges(surface_mesh))
{
halfedge_descriptor hd = halfedge(ed, surface_mesh);
if(is_border(ed, surface_mesh))
{
halfedge_descriptor hd = halfedge(*eb,surface_mesh);
if ( is_border(*eb,surface_mesh) ){
std::cerr << "border" << std::endl;
++nb_sharp_edges;
constraint_hmap[*eb]=true;
constrained_edges[*eb]=std::make_pair(point(source(hd,surface_mesh),surface_mesh),
point(target(hd,surface_mesh),surface_mesh));
constraint_hmap[ed] = true;
constrained_edges[ed] = std::make_pair(point(source(hd, surface_mesh), surface_mesh),
point(target(hd, surface_mesh), surface_mesh));
}
else{
double angle = CGAL::approximate_dihedral_angle(point(target(opposite(hd,surface_mesh),surface_mesh),surface_mesh),
point(target(hd,surface_mesh),surface_mesh),
point(target(next(hd,surface_mesh),surface_mesh),surface_mesh),
point(target(next(opposite(hd,surface_mesh),surface_mesh),surface_mesh),surface_mesh));
if ( CGAL::abs(angle)<100 ){
else
{
double angle = CGAL::approximate_dihedral_angle(point(target(opposite(hd, surface_mesh), surface_mesh), surface_mesh),
point(target(hd, surface_mesh), surface_mesh),
point(target(next(hd, surface_mesh), surface_mesh), surface_mesh),
point(target(next(opposite(hd, surface_mesh), surface_mesh), surface_mesh), surface_mesh));
if(CGAL::abs(angle)<100){
++nb_sharp_edges;
constraint_hmap[*eb]=true;
Point_3 p = point(source(hd,surface_mesh),surface_mesh);
Point_3 q = point(target(hd,surface_mesh),surface_mesh);
constrained_edges[*eb]=std::make_pair(p,q);
constraint_hmap[ed] = true;
Point_3 p = point(source(hd, surface_mesh), surface_mesh);
Point_3 q = point(target(hd, surface_mesh), surface_mesh);
constrained_edges[ed] = std::make_pair(p,q);
cst_output << "2 " << p << " " << q << "\n";
}
}
@ -127,60 +132,60 @@ int main( int argc, char** argv )
// Contract the surface mesh as much as possible
SMS::Count_stop_predicate<Surface_mesh> stop(0);
int r
= SMS::edge_collapse(surface_mesh
,stop
,CGAL::parameters::vertex_index_map(get(CGAL::vertex_external_index, surface_mesh))
int r = SMS::edge_collapse(surface_mesh,
stop,
CGAL::parameters::vertex_index_map(get(CGAL::vertex_external_index, surface_mesh))
.halfedge_index_map(get(CGAL::halfedge_external_index, surface_mesh))
.edge_is_constrained_map(constraints_map)
.get_placement(placement)
);
.get_placement(placement));
std::cout << "\nFinished...\n" << r << " edges removed.\n"
<< num_edges(surface_mesh) << " final edges.\n" ;
std::ofstream os(argc > 2 ? argv[2] : "out.off") ; os << surface_mesh ;
<< num_edges(surface_mesh) << " final edges.\n";
std::ofstream os(argc > 2 ? argv[2] : "out.off"); os << surface_mesh;
std::cout << "Checking sharped edges were preserved...\n";
// check sharp edges were preserved
for(boost::tie(eb,ee) = edges(surface_mesh); eb != ee ; ++eb )
for(edge_descriptor ed : edges(surface_mesh))
{
halfedge_descriptor hd = halfedge(ed, surface_mesh);
if(is_border(ed, surface_mesh))
{
halfedge_descriptor hd = halfedge(*eb,surface_mesh);
if ( is_border(*eb,surface_mesh) ){
--nb_sharp_edges;
assert(
constrained_edges[*eb]==std::make_pair( point(source(hd,surface_mesh),surface_mesh),
point(target(hd,surface_mesh),surface_mesh)));
assert(constrained_edges[ed] == std::make_pair(point(source(hd, surface_mesh), surface_mesh),
point(target(hd, surface_mesh), surface_mesh)));
}
else{
double angle = approximate_dihedral_angle(point(target(opposite(hd,surface_mesh),surface_mesh),surface_mesh),
point(target(hd,surface_mesh),surface_mesh),
point(target(next(hd,surface_mesh),surface_mesh),surface_mesh),
point(target(next(opposite(hd,surface_mesh),surface_mesh),surface_mesh),surface_mesh));
if ( CGAL::abs(angle)<100 ){
else
{
double angle = approximate_dihedral_angle(point(target(opposite(hd, surface_mesh), surface_mesh), surface_mesh),
point(target(hd, surface_mesh), surface_mesh),
point(target(next(hd, surface_mesh), surface_mesh), surface_mesh),
point(target(next(opposite(hd, surface_mesh), surface_mesh), surface_mesh), surface_mesh));
if(CGAL::abs(angle)<100)
{
--nb_sharp_edges;
assert(
constrained_edges[*eb]==std::make_pair( point(source(hd,surface_mesh),surface_mesh),
point(target(hd,surface_mesh),surface_mesh)));
assert(constrained_edges[ed] == std::make_pair(point(source(hd, surface_mesh), surface_mesh),
point(target(hd, surface_mesh), surface_mesh)));
}
}
}
std::cout << "OK\n";
std::cerr << "# sharp edges = " << nb_sharp_edges << std::endl;
std::cout << "Check that no removable edge has been forgotten..." << std::endl;
r = SMS::edge_collapse(surface_mesh
,stop
,CGAL::parameters::vertex_index_map(get(CGAL::vertex_external_index, surface_mesh))
.halfedge_index_map (get(CGAL::halfedge_external_index, surface_mesh))
r = SMS::edge_collapse(surface_mesh,
stop,
CGAL::parameters::vertex_index_map(get(CGAL::vertex_external_index, surface_mesh))
.halfedge_index_map(get(CGAL::halfedge_external_index, surface_mesh))
.edge_is_constrained_map(constraints_map)
.get_placement(placement)
);
assert(r==0);
if ( r==0 )
if(r == 0) {
std::cout << "OK\n";
else{
} else {
std::cout << "ERROR! " << r << " edges removed!\n";
return EXIT_FAILURE;
}

84
Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_constrained_border_polyhedron.cpp
View File

@ -1,7 +1,3 @@
#include <iostream>
#include <fstream>
#include <map>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Polyhedron_3.h>
@ -17,54 +13,56 @@
// Stop-condition policy
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Count_stop_predicate.h>
#include <iostream>
#include <fstream>
#include <map>
typedef CGAL::Simple_cartesian<double> Kernel;
typedef Kernel::Point_3 Point_3;
typedef CGAL::Polyhedron_3<Kernel> Surface_mesh;
namespace SMS = CGAL::Surface_mesh_simplification ;
namespace SMS = CGAL::Surface_mesh_simplification;
//
// BGL property map which indicates whether an edge is marked as non-removable
//
struct Border_is_constrained_edge_map{
struct Border_is_constrained_edge_map
{
const Surface_mesh* sm_ptr;
typedef boost::graph_traits<Surface_mesh>::edge_descriptor key_type;
typedef bool value_type;
typedef value_type reference;
typedef boost::readable_property_map_tag category;
Border_is_constrained_edge_map(const Surface_mesh& sm)
: sm_ptr(&sm)
{}
Border_is_constrained_edge_map(const Surface_mesh& sm) : sm_ptr(&sm) {}
friend bool get(Border_is_constrained_edge_map m, const key_type& edge) {
return CGAL::is_border(edge, *m.sm_ptr);
}
};
//
// Placement class
//
typedef SMS::Constrained_placement<SMS::Midpoint_placement<Surface_mesh>,
Border_is_constrained_edge_map > Placement;
int main( int argc, char** argv )
int main(int argc, char** argv)
{
Surface_mesh surface_mesh;
if (argc!=2){
if(argc != 2)
{
std::cerr << "Usage: " << argv[0] << " input.off\n";
return EXIT_FAILURE;
}
std::ifstream is(argv[1]);
if(!is){
if(!is)
{
std::cerr << "Filename provided is invalid\n";
return EXIT_FAILURE;
}
is >> surface_mesh ;
if (!CGAL::is_triangle_mesh(surface_mesh)){
is >> surface_mesh;
if(!CGAL::is_triangle_mesh(surface_mesh))
{
std::cerr << "Input geometry is not triangulated." << std::endl;
return EXIT_FAILURE;
}
@ -74,56 +72,54 @@ int main( int argc, char** argv )
std::map<Surface_mesh::Halfedge_handle,std::pair<Point_3, Point_3> >constrained_edges;
std::size_t nb_border_edges=0;
for (Surface_mesh::Halfedge_iterator hit=surface_mesh.halfedges_begin(),
for(Surface_mesh::Halfedge_iterator hit=surface_mesh.halfedges_begin(),
hit_end=surface_mesh.halfedges_end();
hit!=hit_end; ++hit )
hit!=hit_end; ++hit)
{
if ( hit->is_border() ){
constrained_edges[hit]=std::make_pair( hit->opposite()->vertex()->point(),
hit->vertex()->point() );
if(hit->is_border())
{
constrained_edges[hit] = std::make_pair(hit->opposite()->vertex()->point(),
hit->vertex()->point());
++nb_border_edges;
}
}
// Contract the surface mesh as much as possible
SMS::Count_stop_predicate<Surface_mesh> stop(0);
Border_is_constrained_edge_map bem(surface_mesh);
// This the actual call to the simplification algorithm.
// The surface mesh and stop conditions are mandatory arguments.
// The index maps are needed because the vertices and edges
// of this surface mesh lack an "id()" field.
int r = SMS::edge_collapse
(surface_mesh
,stop
,CGAL::parameters::vertex_index_map(get(CGAL::vertex_external_index,surface_mesh))
.halfedge_index_map (get(CGAL::halfedge_external_index ,surface_mesh))
int r = SMS::edge_collapse(surface_mesh,
stop,
CGAL::parameters::vertex_index_map(get(CGAL::vertex_external_index,surface_mesh))
.halfedge_index_map(get(CGAL::halfedge_external_index ,surface_mesh))
.edge_is_constrained_map(bem)
.get_placement(Placement(bem))
);
.get_placement(Placement(bem)));
std::cout << "\nFinished...\n" << r << " edges removed.\n"
<< (surface_mesh.size_of_halfedges()/2) << " final edges.\n" ;
<< (surface_mesh.size_of_halfedges()/2) << " final edges.\n";
std::ofstream os( argc > 2 ? argv[2] : "out.off" ) ;
os.precision(17) ;
os << surface_mesh ;
std::ofstream os(argc > 2 ? argv[2] : "out.off");
os.precision(17);
os << surface_mesh;
// now check!
for (Surface_mesh::Halfedge_iterator hit=surface_mesh.halfedges_begin(),
for(Surface_mesh::Halfedge_iterator hit=surface_mesh.halfedges_begin(),
hit_end=surface_mesh.halfedges_end();
hit!=hit_end; ++hit )
hit!=hit_end; ++hit)
{
if(hit->is_border())
{
if (hit->is_border()){
--nb_border_edges;
assert( constrained_edges[hit] ==
std::make_pair( hit->opposite()->vertex()->point(),
hit->vertex()->point() ) );
assert(constrained_edges[hit] == std::make_pair(hit->opposite()->vertex()->point(),
hit->vertex()->point()));
}
}
assert( nb_border_edges==0 );
return EXIT_SUCCESS ;
assert(nb_border_edges==0);
return EXIT_SUCCESS;
}

67
Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_constrained_border_surface_mesh.cpp
View File

@ -1,6 +1,3 @@
#include <iostream>
#include <fstream>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Surface_mesh.h>
@ -16,6 +13,9 @@
// Stop-condition policy
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Count_stop_predicate.h>
#include <iostream>
#include <fstream>
typedef CGAL::Simple_cartesian<double> Kernel;
typedef Kernel::Point_3 Point_3;
@ -25,48 +25,46 @@ typedef boost::graph_traits<Surface_mesh>::edge_descriptor edge_descriptor;
namespace SMS = CGAL::Surface_mesh_simplification;
//
// BGL property map which indicates whether an edge is marked as non-removable
//
struct Border_is_constrained_edge_map{
struct Border_is_constrained_edge_map
{
const Surface_mesh* sm_ptr;
typedef edge_descriptor key_type;
typedef bool value_type;
typedef value_type reference;
typedef boost::readable_property_map_tag category;
Border_is_constrained_edge_map(const Surface_mesh& sm)
: sm_ptr(&sm)
{}
Border_is_constrained_edge_map(const Surface_mesh& sm) : sm_ptr(&sm) {}
friend bool get(Border_is_constrained_edge_map m, const key_type& edge) {
return CGAL::is_border(edge, *m.sm_ptr);
}
};
//
// Placement class
//
typedef SMS::Constrained_placement<SMS::Midpoint_placement<Surface_mesh>,
Border_is_constrained_edge_map > Placement;
int main( int argc, char** argv )
int main(int argc, char** argv)
{
Surface_mesh surface_mesh;
if (argc!=2){
if(argc!=2)
{
std::cerr << "Usage: " << argv[0] << " input.off\n";
return EXIT_FAILURE;
}
std::ifstream is(argv[1]);
if(!is){
if(!is)
{
std::cerr << "Filename provided is invalid\n";
return EXIT_FAILURE;
}
is >> surface_mesh ;
if (!CGAL::is_triangle_mesh(surface_mesh)){
is >> surface_mesh;
if(!CGAL::is_triangle_mesh(surface_mesh))
{
std::cerr << "Input geometry is not triangulated." << std::endl;
return EXIT_FAILURE;
}
@ -76,47 +74,48 @@ int main( int argc, char** argv )
constrained_halfedges = surface_mesh.add_property_map<halfedge_descriptor,std::pair<Point_3, Point_3> >("h:vertices").first;
std::size_t nb_border_edges=0;
for(halfedge_descriptor hd : halfedges(surface_mesh)){
if(CGAL::is_border(hd,surface_mesh)){
constrained_halfedges[hd] = std::make_pair(surface_mesh.point(source(hd,surface_mesh)),
surface_mesh.point(target(hd,surface_mesh)));
for(halfedge_descriptor hd : halfedges(surface_mesh))
{
if(CGAL::is_border(hd, surface_mesh))
{
constrained_halfedges[hd] = std::make_pair(surface_mesh.point(source(hd, surface_mesh)),
surface_mesh.point(target(hd, surface_mesh)));
++nb_border_edges;
}
}
// Contract the surface mesh as much as possible
SMS::Count_stop_predicate<Surface_mesh> stop(0);
Border_is_constrained_edge_map bem(surface_mesh);
// This the actual call to the simplification algorithm.
// The surface mesh and stop conditions are mandatory arguments.
int r = SMS::edge_collapse
(surface_mesh
,stop
,CGAL::parameters::edge_is_constrained_map(bem)
.get_placement(Placement(bem))
);
int r = SMS::edge_collapse(surface_mesh,
stop,
CGAL::parameters::edge_is_constrained_map(bem)
.get_placement(Placement(bem)));
std::cout << "\nFinished...\n" << r << " edges removed.\n"
<< surface_mesh.number_of_edges() << " final edges.\n";
std::ofstream os( argc > 2 ? argv[2] : "out.off" );
std::ofstream os(argc > 2 ? argv[2] : "out.off");
os.precision(17);
os << surface_mesh;
// now check!
for(halfedge_descriptor hd : halfedges(surface_mesh)){
if(CGAL::is_border(hd,surface_mesh)){
for(halfedge_descriptor hd : halfedges(surface_mesh))
{
if(CGAL::is_border(hd,surface_mesh))
{
--nb_border_edges;
if(constrained_halfedges[hd] != std::make_pair(surface_mesh.point(source(hd,surface_mesh)),
surface_mesh.point(target(hd,surface_mesh)))){
if(constrained_halfedges[hd] != std::make_pair(surface_mesh.point(source(hd, surface_mesh)),
surface_mesh.point(target(hd, surface_mesh))))
{
std::cerr << "oops. send us a bug report\n";
}
}
}
assert( nb_border_edges==0 );
assert(nb_border_edges==0);
return EXIT_SUCCESS;
}

6
Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_enriched_polyhedron.cpp
View File

@ -22,13 +22,13 @@ typedef boost::graph_traits<Surface_mesh>::vertex_descriptor vertex_descriptor;
namespace SMS = CGAL::Surface_mesh_simplification;
int main( int argc, char** argv )
int main(int argc, char** argv)
{
Surface_mesh surface_mesh;
std::ifstream is(argv[1]);
is >> surface_mesh;
if (!CGAL::is_triangle_mesh(surface_mesh)){
if(!CGAL::is_triangle_mesh(surface_mesh)){
std::cerr << "Input geometry is not triangulated." << std::endl;
return EXIT_FAILURE;
}
@ -64,7 +64,7 @@ int main( int argc, char** argv )
std::cout << "\nFinished...\n" << r << " edges removed.\n"
<< (surface_mesh.size_of_halfedges()/2) << " final edges.\n";
std::ofstream os( argc > 2 ? argv[2] : "out.off" );
std::ofstream os(argc > 2 ? argv[2] : "out.off");
os.precision(17);
os << surface_mesh;

8
Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_linear_cell_complex.cpp
View File

@ -12,11 +12,11 @@ typedef CGAL::Simple_cartesian<double> Kernel;
typedef CGAL::Linear_cell_complex_traits<3, Kernel> MyTraits;
typedef CGAL::Linear_cell_complex_for_bgl_combinatorial_map_helper
<2, 3, MyTraits>::type LCC;
namespace SMS = CGAL::Surface_mesh_simplification ;
namespace SMS = CGAL::Surface_mesh_simplification;
int main( int argc, char** argv )
int main(int argc, char** argv)
{
if (argc<2 || argc>3)
if(argc<2 || argc>3)
{
std::cout<<"Usage: simplification_Linear_cell_complex inofffile [outofffile]"<<std::endl;
return EXIT_FAILURE;
@ -45,7 +45,7 @@ int main( int argc, char** argv )
);
std::cout << "\nFinished...\n" << r << " edges removed.\n"
<< (lcc.number_of_darts()/2) << " final edges.\n" ;
<< (lcc.number_of_darts()/2) << " final edges.\n";
lcc.display_characteristics(std::cout)<<", is_valid="<<CGAL::is_valid(lcc)<<std::endl;

6
Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_polyhedron.cpp
View File

@ -16,13 +16,13 @@ typedef CGAL::Polyhedron_3<Kernel> Surface_mesh;
namespace SMS = CGAL::Surface_mesh_simplification;
int main( int argc, char** argv )
int main(int argc, char** argv)
{
Surface_mesh surface_mesh;
std::ifstream is(argv[1]);
is >> surface_mesh;
if (!CGAL::is_triangle_mesh(surface_mesh)){
if(!CGAL::is_triangle_mesh(surface_mesh)){
std::cerr << "Input geometry is not triangulated." << std::endl;
return EXIT_FAILURE;
}
@ -46,7 +46,7 @@ int main( int argc, char** argv )
std::cout << "\nFinished...\n" << r << " edges removed.\n"
<< (surface_mesh.size_of_halfedges()/2) << " final edges.\n";
std::ofstream os( argc > 2 ? argv[2] : "out.off" );
std::ofstream os(argc > 2 ? argv[2] : "out.off");
os.precision(17);
os << surface_mesh;

6
Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_surface_mesh.cpp
View File

@ -13,13 +13,13 @@ typedef CGAL::Surface_mesh<Point_3> Surface_mesh;
namespace SMS = CGAL::Surface_mesh_simplification;
int main( int argc, char** argv )
int main(int argc, char** argv)
{
Surface_mesh surface_mesh;
std::ifstream is(argv[1]);
is >> surface_mesh;
if (!CGAL::is_triangle_mesh(surface_mesh)){
if(!CGAL::is_triangle_mesh(surface_mesh)){
std::cerr << "Input geometry is not triangulated." << std::endl;
return EXIT_FAILURE;
}
@ -33,7 +33,7 @@ int main( int argc, char** argv )
std::cout << "\nFinished...\n" << r << " edges removed.\n"
<< surface_mesh.number_of_edges() << " final edges.\n";
std::ofstream os( argc > 2 ? argv[2] : "out.off" );
std::ofstream os(argc > 2 ? argv[2] : "out.off");
os.precision(17);
os << surface_mesh;

26
Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_visitor_surface_mesh.cpp
View File

@ -45,16 +45,16 @@ struct Stats
std::size_t processed;
std::size_t collapsed;
std::size_t non_collapsable;
std::size_t cost_uncomputable ;
std::size_t cost_uncomputable;
std::size_t placement_uncomputable;
};
struct My_visitor : SMS::Edge_collapse_visitor_base<Surface_mesh>
{
My_visitor( Stats* s) : stats(s){}
My_visitor(Stats* s) : stats(s){}
// Called during the collecting phase for each edge collected.
void OnCollected( Profile const&, boost::optional<double> const& )
void OnCollected(Profile const&, boost::optional<double> const&)
{
++ stats->collected;
std::cerr << "\rEdges collected: " << stats->collected << std::flush;
@ -62,41 +62,41 @@ struct My_visitor : SMS::Edge_collapse_visitor_base<Surface_mesh>
// Called during the processing phase for each edge selected.
// If cost is absent the edge won't be collapsed.
void OnSelected(Profile const&
void OnSelected(const Profile&
,boost::optional<double> cost
,std::size_t initial
,std::size_t current
)
{
++ stats->processed;
if ( !cost )
if(!cost)
++ stats->cost_uncomputable;
if ( current == initial )
if(current == initial)
std::cerr << "\n" << std::flush;
std::cerr << "\r" << current << std::flush;
}
// Called during the processing phase for each edge being collapsed.
// If placement is absent the edge is left uncollapsed.
void OnCollapsing(Profile const&
void OnCollapsing(const Profile&
,boost::optional<Point> placement
)
{
if ( !placement )
if(!placement)
++ stats->placement_uncomputable;
}
// Called for each edge which failed the so called link-condition,
// that is, which cannot be collapsed because doing so would
// turn the surface mesh into a non-manifold.
void OnNonCollapsable( Profile const& )
void OnNonCollapsable(Profile const&)
{
++ stats->non_collapsable;
}
// Called after each edge has been collapsed
void OnCollapsed( Profile const&, vertex_descriptor )
void OnCollapsed(Profile const&, vertex_descriptor)
{
++ stats->collapsed;
}
@ -105,13 +105,13 @@ struct My_visitor : SMS::Edge_collapse_visitor_base<Surface_mesh>
};
int main( int argc, char** argv )
int main(int argc, char** argv)
{
Surface_mesh surface_mesh;
std::ifstream is(argv[1]);
is >> surface_mesh;
if (!CGAL::is_triangle_mesh(surface_mesh)){
if(!CGAL::is_triangle_mesh(surface_mesh)){
std::cerr << "Input geometry is not triangulated." << std::endl;
return EXIT_FAILURE;
}
@ -147,7 +147,7 @@ int main( int argc, char** argv )
std::cout << "\nFinished...\n" << r << " edges removed.\n"
<< surface_mesh.number_of_edges() << " final edges.\n";
std::ofstream os( argc > 2 ? argv[2] : "out.off" );
std::ofstream os(argc > 2 ? argv[2] : "out.off");
os.precision(17);
os << surface_mesh;

198
Surface_mesh_simplification/include/CGAL/Cartesian/MatrixC33.h
View File

@ -22,8 +22,6 @@
#include <CGAL/license/Surface_mesh_simplification.h>
#include <CGAL/Vector_3.h>
#include <CGAL/determinant.h>
#include <CGAL/number_utils.h>
#include <CGAL/Null_matrix.h>
@ -36,124 +34,116 @@ template <class R_>
class MatrixC33
{
public:
typedef R_ R;
typedef R_ R ;
typedef typename R::FT RT ;
typedef typename R::FT RT;
typedef typename R::Vector_3 Vector_3;
MatrixC33 ( Null_matrix )
:
mR0(NULL_VECTOR)
,mR1(NULL_VECTOR)
,mR2(NULL_VECTOR)
MatrixC33(Null_matrix)
: mR0(NULL_VECTOR),
mR1(NULL_VECTOR),
mR2(NULL_VECTOR)
{}
MatrixC33 ( RT const& r0x, RT const& r0y, RT const& r0z
, RT const& r1x, RT const& r1y, RT const& r1z
, RT const& r2x, RT const& r2y, RT const& r2z
)
:
mR0(r0x,r0y,r0z)
,mR1(r1x,r1y,r1z)
,mR2(r2x,r2y,r2z)
MatrixC33(const RT& r0x, const RT& r0y, const RT& r0z,
const RT& r1x, const RT& r1y, const RT& r1z,
const RT& r2x, const RT& r2y, const RT& r2z)
: mR0(r0x,r0y,r0z),
mR1(r1x,r1y,r1z),
mR2(r2x,r2y,r2z)
{}
MatrixC33 ( Vector_3 const& r0, Vector_3 const& r1, Vector_3 const& r2 )
:
mR0(r0)
,mR1(r1)
,mR2(r2)
MatrixC33(const Vector_3& r0, const Vector_3& r1, const Vector_3& r2)
: mR0(r0),
mR1(r1),
mR2(r2)
{}
Vector_3 const& r0() const { return mR0; }
Vector_3 const& r1() const { return mR1; }
Vector_3 const& r2() const { return mR2; }
const Vector_3& r0() const { return mR0; }
const Vector_3& r1() const { return mR1; }
const Vector_3& r2() const { return mR2; }
Vector_3& r0() { return mR0; }
Vector_3& r1() { return mR1; }
Vector_3& r2() { return mR2; }
Vector_3 const& operator[] ( int row ) const { return row == 0 ? mR0 : ( row == 1 ? mR1 : mR2 ) ; }
Vector_3& operator[] ( int row ) { return row == 0 ? mR0 : ( row == 1 ? mR1 : mR2 ) ; }
const Vector_3& operator[](int row) const { return row == 0 ? mR0 : (row == 1 ? mR1 : mR2); }
Vector_3& operator[](int row) { return row == 0 ? mR0 : (row == 1 ? mR1 : mR2); }
MatrixC33& operator+= ( MatrixC33 const& m )
MatrixC33& operator+=(const MatrixC33& m)
{
mR0 = mR0 + m.r0() ;
mR1 = mR1 + m.r1() ;
mR2 = mR2 + m.r2() ;
return *this ;
mR0 = mR0 + m.r0();
mR1 = mR1 + m.r1();
mR2 = mR2 + m.r2();
return *this;
}
MatrixC33& operator-= ( MatrixC33 const& m )
MatrixC33& operator-=(const MatrixC33& m)
{
mR0 = mR0 - m.r0() ;
mR1 = mR1 - m.r1() ;
mR2 = mR2 - m.r2() ;
return *this ;
mR0 = mR0 - m.r0();
mR1 = mR1 - m.r1();
mR2 = mR2 - m.r2();
return *this;
}
MatrixC33& operator*= ( RT const& c )
MatrixC33& operator*=(const RT& c)
{
mR0 = mR0 * c ;
mR1 = mR1 * c ;
mR2 = mR2 * c ;
return *this ;
mR0 = mR0 * c;
mR1 = mR1 * c;
mR2 = mR2 * c;
return *this;
}
MatrixC33& operator/= ( RT const& c )
MatrixC33& operator/=(const RT& c)
{
mR0 = mR0 / c ;
mR1 = mR1 / c ;
mR2 = mR2 / c ;
return *this ;
mR0 = mR0 / c;
mR1 = mR1 / c;
mR2 = mR2 / c;
return *this;
}
friend MatrixC33 operator+ ( MatrixC33 const& a, MatrixC33 const& b )
friend MatrixC33 operator+(const MatrixC33& a, const MatrixC33& b)
{
return MatrixC33(a.r0()+b.r0()
,a.r1()+b.r1()
,a.r2()+b.r2()
);
return MatrixC33(a.r0() + b.r0(),
a.r1() + b.r1(),
a.r2() + b.r2());
}
friend MatrixC33 operator- ( MatrixC33 const& a, MatrixC33 const& b )
friend MatrixC33 operator-(const MatrixC33& a, const MatrixC33& b)
{
return MatrixC33(a.r0()-b.r0()
,a.r1()-b.r1()
,a.r2()-b.r2()
);
return MatrixC33(a.r0() - b.r0(),
a.r1() - b.r1(),
a.r2() - b.r2());
}
friend MatrixC33 operator* ( MatrixC33 const& m, RT const& c )
friend MatrixC33 operator*(const MatrixC33& m, const RT& c)
{
return MatrixC33(m.r0()*c,m.r1()*c,m.r2()*c);
return MatrixC33(m.r0()*c, m.r1()*c, m.r2()*c);
}
friend MatrixC33 operator* ( RT const& c, MatrixC33 const& m )
friend MatrixC33 operator*(const RT& c, const MatrixC33& m)
{
return MatrixC33(m.r0()*c,m.r1()*c,m.r2()*c);
return MatrixC33(m.r0()*c, m.r1()*c, m.r2()*c);
}
friend MatrixC33 operator/ ( MatrixC33 const& m, RT const& c )
friend MatrixC33 operator/(const MatrixC33& m, const RT& c)
{
return MatrixC33(m.r0()/c,m.r1()/c,m.r2()/c);
return MatrixC33(m.r0()/c, m.r1()/c, m.r2()/c);
}
friend Vector_3 operator* ( MatrixC33 const& m, Vector_3 const& v )
friend Vector_3 operator*(const MatrixC33& m, const Vector_3& v)
{
return Vector_3(m.r0()*v,m.r1()*v,m.r2()*v);
return Vector_3(m.r0()*v, m.r1()*v, m.r2()*v);
}
friend Vector_3 operator* ( Vector_3 const& v, MatrixC33 const& m )
friend Vector_3 operator*(const Vector_3& v, const MatrixC33& m)
{
return Vector_3(v*m.r0(),v*m.r1(),v*m.r2());
return Vector_3(v*m.r0(), v*m.r1(), v*m.r2());
}
RT determinant() const
{
return CGAL::determinant(r0().x(),r0().y(),r0().z()
,r1().x(),r1().y(),r1().z()
,r2().x(),r2().y(),r2().z()
);
return CGAL::determinant(r0().x(), r0().y(), r0().z(),
r1().x(), r1().y(), r1().z(),
r2().x(), r2().y(), r2().z());
}
MatrixC33& transpose()
@ -161,38 +151,37 @@ public:
mR0 = Vector_3(r0().x(),r1().x(),r2().x());
mR1 = Vector_3(r0().y(),r1().y(),r2().y());
mR2 = Vector_3(r0().z(),r1().z(),r2().z());
return *this ;
return *this;
}
private:
Vector_3 mR0 ;
Vector_3 mR1 ;
Vector_3 mR2 ;
} ;
Vector_3 mR0;
Vector_3 mR1;
Vector_3 mR2;
};
template<class R>
inline
MatrixC33<R> direct_product ( Vector_3<R> const& u, Vector_3<R> const& v )
MatrixC33<R> direct_product(Vector_3<R> const& u, Vector_3<R> const& v)
{
return MatrixC33<R>( v * u.x()
, v * u.y()
, v * u.z()
) ;
return MatrixC33<R>(v * u.x(),
v * u.y(),
v * u.z());
}
template<class R>
MatrixC33<R> transposed_matrix ( MatrixC33<R> const& m )
MatrixC33<R> transposed_matrix(const MatrixC33<R>& m)
{
MatrixC33<R> copy = m ;
MatrixC33<R> copy = m;
copy.Transpose();
return copy ;
return copy;
}
template<class R>
MatrixC33<R> cofactors_matrix ( MatrixC33<R> const& m )
MatrixC33<R> cofactors_matrix(const MatrixC33<R>& m)
{
typedef typename R::RT RT ;
typedef typename R::RT RT;
RT c00 = determinant(m.r1().y(),m.r1().z(),m.r2().y(),m.r2().z());
RT c01 = -determinant(m.r1().x(),m.r1().z(),m.r2().x(),m.r2().z());
@ -206,32 +195,29 @@ MatrixC33<R> cofactors_matrix ( MatrixC33<R> const& m )
RT c21 = -determinant(m.r0().x(),m.r0().z(),m.r1().x(),m.r1().z());
RT c22 = determinant(m.r0().x(),m.r0().y(),m.r1().x(),m.r1().y());
return MatrixC33<R>(c00,c01,c02
,c10,c11,c12
,c20,c21,c22
);
return MatrixC33<R>(c00,c01,c02,
c10,c11,c12,
c20,c21,c22);
}
template<class R>
MatrixC33<R> adjoint_matrix ( MatrixC33<R> const& m )
MatrixC33<R> adjoint_matrix(const MatrixC33<R>& m)
{
return cofactors_matrix(m).transpose() ;
return cofactors_matrix(m).transpose();
}
template<class R>
boost::optional< MatrixC33<R> > inverse_matrix ( MatrixC33<R> const& m )
boost::optional< MatrixC33<R> > inverse_matrix(const MatrixC33<R>& m)
{
typedef typename R::RT RT ;
typedef typename R::RT RT;
typedef MatrixC33<R> Matrix;
typedef boost::optional<Matrix> result_type;
typedef MatrixC33<R> Matrix ;
typedef boost::optional<Matrix> result_type ;
result_type rInverse ;
result_type rInverse;
RT det = m.determinant();
if ( ! CGAL_NTS is_zero(det) )
if(! CGAL_NTS is_zero(det))
{
RT c00 = (m.r1().y()*m.r2().z() - m.r1().z()*m.r2().y()) / det;
RT c01 = (m.r2().y()*m.r0().z() - m.r0().y()*m.r2().z()) / det;
@ -245,14 +231,12 @@ boost::optional< MatrixC33<R> > inverse_matrix ( MatrixC33<R> const& m )
RT c21 = (m.r2().x()*m.r0().y() - m.r0().x()*m.r2().y()) / det;
RT c22 = (m.r0().x()*m.r1().y() - m.r0().y()*m.r1().x()) / det;
rInverse = result_type( Matrix(c00,c01,c02
,c10,c11,c12
,c20,c21,c22
)
) ;
rInverse = result_type(Matrix(c00,c01,c02,
c10,c11,c12,
c20,c21,c22));
}
return rInverse ;
return rInverse;
}
} //namespace CGAL

11
Surface_mesh_simplification/include/CGAL/Null_matrix.h
View File

@ -22,19 +22,12 @@
#include <CGAL/license/Surface_mesh_simplification.h>
namespace CGAL {
class Null_matrix {};
static const Null_matrix NULL_MATRIX = Null_matrix() ;
static const Null_matrix NULL_MATRIX = Null_matrix();
} //namespace CGAL
#endif // CGAL_NULL_MATRIX_H //
// EOF //
#endif // CGAL_NULL_MATRIX_H

153
Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Detail/Common.h
View File

@ -18,16 +18,15 @@
// Author(s) : Fernando Cacciola <fernando.cacciola@geometryfactory.com>
//
#ifndef CGAL_SURFACE_MESH_SIMPLIFICATION_DETAIL_COMMON_H
#define CGAL_SURFACE_MESH_SIMPLIFICATION_DETAIL_COMMON_H 1
#define CGAL_SURFACE_MESH_SIMPLIFICATION_DETAIL_COMMON_H
#include <CGAL/license/Surface_mesh_simplification.h>
#include <functional>
#include <utility>
#include <vector>
#include <vector>
#include <set>
#include <CGAL/algorithm.h>
#include <CGAL/Cartesian/MatrixC33.h>
#include <CGAL/Modifiable_priority_queue.h>
#include <CGAL/boost/graph/properties.h>
#include <CGAL/boost/graph/iterator.h>
#include <boost/config.hpp>
#include <boost/shared_ptr.hpp>
@ -42,91 +41,84 @@
#include <boost/graph/properties.hpp>
#include <boost/graph/adjacency_list.hpp>
#include <CGAL/algorithm.h>
#include <CGAL/Cartesian/MatrixC33.h>
#include <CGAL/Modifiable_priority_queue.h>
#include <CGAL/boost/graph/properties.h>
#include <CGAL/boost/graph/iterator.h>
#include <functional>
#include <utility>
#include <vector>
#include <vector>
#include <set>
namespace CGAL {
namespace Surface_mesh_simplification {
namespace Surface_mesh_simplification
{
using boost::num_edges;
using boost::num_vertices;
using boost::edges;
using boost::out_edges;
using boost::in_edges;
using boost::source;
using boost::target;
using boost::num_edges ;
using boost::num_vertices ;
using boost::edges ;
using boost::out_edges ;
using boost::in_edges ;
using boost::source ;
using boost::target ;
using boost::shared_ptr;
using boost::optional;
using boost::none;
using boost::put_get_helper;
using boost::get;
using boost::put;
using boost::addressof;
using boost::shared_ptr ;
using boost::optional ;
using boost::none ;
using boost::put_get_helper ;
using boost::get ;
using boost::put ;
using boost::addressof ;
using namespace boost::tuples ;
using namespace boost::tuples;
template<class Handle>
inline bool handle_assigned( Handle h ) { Handle null ; return h != null ; }
inline bool handle_assigned(Handle h) { Handle null; return h != null; }
template<class Iterator, class Handle>
bool handle_exists ( Iterator begin, Iterator end, Handle h )
bool handle_exists(Iterator begin, Iterator end, Handle h)
{
if ( handle_assigned(h) )
if(handle_assigned(h))
{
while ( begin != end )
if ( begin++ == h )
return true ;
while(begin != end)
if(begin++ == h)
return true;
}
return false ;
return false;
}
template <class TM>
struct No_constrained_edge_map{
struct No_constrained_edge_map
{
typedef typename boost::graph_traits<TM>::edge_descriptor key_type;
typedef bool value_type;
typedef value_type reference;
typedef boost::readable_property_map_tag category;
friend bool get(No_constrained_edge_map, key_type) {
return false;
}
friend bool get(No_constrained_edge_map, key_type) { return false; }
};
} // namespace Surface_mesh_simplification
template<class N>
inline std::string n_to_string( N const& n )
{
return boost::str( boost::format("%|5.19g|") % n ) ;
inline std::string n_to_string(const N& n) {
return boost::str(boost::format("%|5.19g|") % n);
}
template<class XYZ>
inline std::string xyz_to_string( XYZ const& xyz )
{
return boost::str( boost::format("(%|5.19g|,%|5.19g|,%|5.19g|)") % xyz.x() % xyz.y() % xyz.z() ) ;
inline std::string xyz_to_string(const XYZ& xyz) {
return boost::str(boost::format("(%|5.19g|,%|5.19g|,%|5.19g|)") % xyz.x() % xyz.y() % xyz.z());
}
template<class Matrix>
inline std::string matrix_to_string( Matrix const& m )
{
return boost::str( boost::format("[%1%|%2%|%3%]") % xyz_to_string(m.r0()) % xyz_to_string(m.r1()) % xyz_to_string(m.r2()) ) ;
inline std::string matrix_to_string(const Matrix& m) {
return boost::str(boost::format("[%1%|%2%|%3%]") % xyz_to_string(m.r0()) % xyz_to_string(m.r1()) % xyz_to_string(m.r2()));
}
template<class T>
inline std::string optional_to_string( boost::optional<T> const& o )
{
if ( o )
return boost::str( boost::format("%1%") % *o ) ;
inline std::string optional_to_string(const boost::optional<T>& o) {
if(o)
return boost::str(boost::format("%1%") % *o);
else return std::string("NONE");
}
} //namespace CGAL
} // namespace CGAL
#if defined(CGAL_SURFACE_SIMPLIFICATION_ENABLE_TRACE) \
|| defined(CGAL_SURFACE_SIMPLIFICATION_ENABLE_LT_TRACE)
@ -134,52 +126,45 @@ inline std::string optional_to_string( boost::optional<T> const& o )
#endif
#ifdef CGAL_SMS_ENABLE_TRACE
#include<string>
#include<iostream>
#include<sstream>
# include<string>
# include<iostream>
# include<sstream>
namespace internal { namespace { bool cgal_enable_sms_trace = true ; } }
# define CGAL_SMS_TRACE_IMPL(m) \
if ( ::internal::cgal_enable_sms_trace ) { \
std::ostringstream ss ; ss << m ; std::string s = ss.str(); \
namespace internal { namespace { bool cgal_enable_sms_trace = true; } }
#define CGAL_SMS_TRACE_IMPL(m) \
if(::internal::cgal_enable_sms_trace) { \
std::ostringstream ss; ss << m; std::string s = ss.str(); \
/*Surface_simplification_external_trace(s)*/ std::cerr << s << std::endl; \
}
# define CGAL_SMS_DEBUG_CODE(code) code
#define CGAL_SMS_DEBUG_CODE(code) code
#else
# define CGAL_SMS_DEBUG_CODE(code)
#define CGAL_SMS_DEBUG_CODE(code)
#endif
#ifdef CGAL_SURFACE_SIMPLIFICATION_ENABLE_LT_TRACE
# define CGAL_SMS_LT_TRACE(l,m) if ( (l) <= CGAL_SURFACE_SIMPLIFICATION_ENABLE_LT_TRACE ) CGAL_SMS_TRACE_IMPL(m)
#define CGAL_SMS_LT_TRACE(l,m) if((l) <= CGAL_SURFACE_SIMPLIFICATION_ENABLE_LT_TRACE) CGAL_SMS_TRACE_IMPL(m)
#else
# define CGAL_SMS_LT_TRACE(l,m)
#define CGAL_SMS_LT_TRACE(l,m)
#endif
#ifdef CGAL_SURFACE_SIMPLIFICATION_ENABLE_TRACE
# define CGAL_SMS_TRACE_IF(c,l,m) if ( (c) && ( (l) <= CGAL_SURFACE_SIMPLIFICATION_ENABLE_TRACE) ) CGAL_SMS_TRACE_IMPL(m)
# define CGAL_SMS_TRACE(l,m) if ( (l) <= CGAL_SURFACE_SIMPLIFICATION_ENABLE_TRACE ) CGAL_SMS_TRACE_IMPL(m)
#define CGAL_SMS_TRACE_IF(c,l,m) if((c) && ((l) <= CGAL_SURFACE_SIMPLIFICATION_ENABLE_TRACE)) CGAL_SMS_TRACE_IMPL(m)
#define CGAL_SMS_TRACE(l,m) if((l) <= CGAL_SURFACE_SIMPLIFICATION_ENABLE_TRACE) CGAL_SMS_TRACE_IMPL(m)
#else
# define CGAL_SMS_TRACE_IF(c,l,m)
# define CGAL_SMS_TRACE(l,m)
#define CGAL_SMS_TRACE_IF(c,l,m)
#define CGAL_SMS_TRACE(l,m)
#endif
#undef CGAL_SMS_ENABLE_TRACE
#ifdef CGAL_TESTING_SURFACE_MESH_SIMPLIFICATION
# define CGAL_SURF_SIMPL_TEST_assertion(EX) CGAL_assertion(EX)
# define CGAL_SURF_SIMPL_TEST_assertion_msg(EX,MSG) CGAL_assertion_msg(EX,MSG)
# define CGAL_SURF_SIMPL_TEST_assertion_code(CODE) CGAL_assertion_code(CODE)
#define CGAL_SURF_SIMPL_TEST_assertion(EX) CGAL_assertion(EX)
#define CGAL_SURF_SIMPL_TEST_assertion_msg(EX,MSG) CGAL_assertion_msg(EX,MSG)
#define CGAL_SURF_SIMPL_TEST_assertion_code(CODE) CGAL_assertion_code(CODE)
#else
# define CGAL_SURF_SIMPL_TEST_assertion(EX)
# define CGAL_SURF_SIMPL_TEST_assertion_msg(EX,MSG)
# define CGAL_SURF_SIMPL_TEST_assertion_code(CODE)
#define CGAL_SURF_SIMPL_TEST_assertion(EX)
#define CGAL_SURF_SIMPL_TEST_assertion_msg(EX,MSG)
#define CGAL_SURF_SIMPL_TEST_assertion_code(CODE)
#endif
#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_DETAIL_COMMON_H //
// EOF //

1285
Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Detail/Edge_collapse.h
View File

File diff suppressed because it is too large Load Diff

939
Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Detail/Edge_collapse_impl.h
View File

@ -1,939 +0,0 @@
// Copyright (c) 2006 GeometryFactory (France). All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
// You can redistribute it and/or modify it under the terms of the GNU
// General Public License as published by the Free Software Foundation,
// either version 3 of the License, or (at your option) any later version.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0+
//
// Author(s) : Fernando Cacciola <fernando.cacciola@geometryfactory.com>
//
#ifndef CGAL_SURFACE_MESH_SIMPLIFICATION_DETAIL_EDGE_COLLAPSE_IMPL_H
#define CGAL_SURFACE_MESH_SIMPLIFICATION_DETAIL_EDGE_COLLAPSE_IMPL_H
#include <CGAL/license/Surface_mesh_simplification.h>
namespace CGAL {
namespace Surface_mesh_simplification
{
template<class M, class SP, class VIM, class VPM,class EIM, class ECTM, class CF, class PF, class V>
EdgeCollapse<M,SP,VIM,VPM,EIM,ECTM,CF,PF,V>::EdgeCollapse( TM& aSurface
, ShouldStop const& aShould_stop
, VertexIndexMap const& aVertex_index_map
, VertexPointMap const& aVertex_point_map
, EdgeIndexMap const& aEdge_index_map
, EdgeIsConstrainedMap const& aEdge_is_constrained_map
, GetCost const& aGet_cost
, GetPlacement const& aGet_placement
, VisitorT aVisitor
)
:
mSurface (aSurface)
,Should_stop (aShould_stop)
,Vertex_index_map (aVertex_index_map)
,Vertex_point_map (aVertex_point_map)
,Edge_index_map (aEdge_index_map)
,Edge_is_constrained_map (aEdge_is_constrained_map)
,Get_cost (aGet_cost)
,Get_placement (aGet_placement)
,Visitor (aVisitor)
,m_has_border (false)
{
const FT cMaxDihedralAngleCos = std::cos( 1.0 * CGAL_PI / 180.0 ) ;
mcMaxDihedralAngleCos2 = cMaxDihedralAngleCos * cMaxDihedralAngleCos ;
halfedge_iterator eb, ee ;
for ( boost::tie(eb,ee) = halfedges(mSurface); eb!=ee; ++eb )
{
halfedge_descriptor ed = *eb;
if(is_border(ed)){
m_has_border = true;
break;
}
}
CGAL_SMS_TRACE(0,"EdgeCollapse of TM with " << (num_edges(aSurface)/2) << " edges" );
CGAL_SMS_DEBUG_CODE ( mStep = 0 ; )
#ifdef CGAL_SURFACE_SIMPLIFICATION_ENABLE_TRACE
vertex_iterator vb, ve ;
for ( boost::tie(vb,ve) = vertices(mSurface) ; vb != ve ; ++ vb )
CGAL_SMS_TRACE(1, vertex_to_string(*vb) ) ;
for ( boost::tie(eb,ee) = halfedges(mSurface); eb!=ee; ++eb )
CGAL_SMS_TRACE(1, edge_to_string(*eb) ) ;
#endif
}
template<class M,class SP, class VIM, class VPM,class EIM,class ECTM, class CF,class PF,class V>
int EdgeCollapse<M,SP,VIM,VPM,EIM,ECTM,CF,PF,V>::run()
{
CGAL_SURF_SIMPL_TEST_assertion( mSurface.is_valid() && mSurface.is_pure_triangle() ) ;
Visitor.OnStarted(mSurface);
// First collect all candidate edges in a PQ
Collect();
// Then proceed to collapse each edge in turn
Loop();
CGAL_SMS_TRACE(0,"Finished: " << (mInitialEdgeCount - mCurrentEdgeCount) << " edges removed." ) ;
int r = (int)(mInitialEdgeCount - mCurrentEdgeCount) ;
Visitor.OnFinished(mSurface);
return r ;
}
template<class M,class SP, class VIM, class VPM,class EIM,class ECTM, class CF,class PF,class V>
void EdgeCollapse<M,SP,VIM,VPM,EIM,ECTM,CF,PF,V>::Collect()
{
CGAL_SMS_TRACE(0,"Collecting edges...");
//
// Loop over all the _undirected_ edges in the surface putting them in the PQ
//
Equal_3 equal_points = Traits().equal_3_object();
size_type lSize = num_edges(mSurface);
mInitialEdgeCount = mCurrentEdgeCount = static_cast<size_type>(
std::distance( boost::begin(edges(mSurface)),
boost::end(edges(mSurface)) ) );;
mEdgeDataArray.reset( new Edge_data[lSize] ) ;
mPQ.reset( new PQ (lSize, Compare_cost(this), edge_id(this) ) ) ;
std::size_t id = 0 ;
CGAL_SURF_SIMPL_TEST_assertion_code ( size_type lInserted = 0 ) ;
CGAL_SURF_SIMPL_TEST_assertion_code ( size_type lNotInserted = 0 ) ;
std::set<halfedge_descriptor> zero_length_edges;
edge_iterator eb, ee ;
for ( boost::tie(eb,ee) = edges(mSurface); eb!=ee; ++eb, id+=2 )
{
halfedge_descriptor lEdge = halfedge(*eb,mSurface) ;
if ( is_constrained(lEdge) ) continue;//no not insert constrainted edges
// AF CGAL_assertion( get_halfedge_id(lEdge) == id ) ;
// AF CGAL_assertion( get_halfedge_id(opposite(lEdge,mSurface)) == id+1 ) ;
Profile const& lProfile = create_profile(lEdge);
if ( !equal_points(lProfile.p0(),lProfile.p1()) )
{
Edge_data& lData = get_data(lEdge);
lData.cost() = get_cost(lProfile) ;
insert_in_PQ(lEdge,lData);
Visitor.OnCollected(lProfile,lData.cost());
CGAL_SURF_SIMPL_TEST_assertion_code ( ++ lInserted ) ;
}
else
{
zero_length_edges.insert(primary_edge(lEdge));
CGAL_SURF_SIMPL_TEST_assertion_code ( ++ lNotInserted ) ;
}
CGAL_SMS_TRACE(2,edge_to_string(lEdge));
}
CGAL_SURF_SIMPL_TEST_assertion ( lInserted + lNotInserted == mInitialEdgeCount ) ;
for (typename std::set<halfedge_descriptor>::iterator it=zero_length_edges.begin(),
it_end=zero_length_edges.end();it!=it_end;++it)
{
Profile const& lProfile = create_profile(*it);
if (!Is_collapse_topologically_valid(lProfile) ) continue;
// edges of length 0 removed no longer need to be treated
if ( lProfile.left_face_exists() )
{
halfedge_descriptor lEdge_to_remove = is_constrained(lProfile.vL_v0()) ?
primary_edge(lProfile.v1_vL()) :
primary_edge(lProfile.vL_v0()) ;
zero_length_edges.erase( lEdge_to_remove );
Edge_data& lData = get_data(lEdge_to_remove) ;
if ( lData.is_in_PQ() ){
CGAL_SMS_TRACE(2,"Removing E" << get(Edge_index_map,lEdge_to_remove) << " from PQ" );
remove_from_PQ(lEdge_to_remove,lData);
}
--mCurrentEdgeCount;
}
if ( lProfile.right_face_exists() )
{
halfedge_descriptor lEdge_to_remove = is_constrained(lProfile.vR_v1()) ?
primary_edge(lProfile.v0_vR()) :
primary_edge(lProfile.vR_v1()) ;
zero_length_edges.erase( lEdge_to_remove );
Edge_data& lData = get_data(lEdge_to_remove) ;
if ( lData.is_in_PQ() ){
CGAL_SMS_TRACE(2,"Removing E" << get(Edge_index_map,lEdge_to_remove) << " from PQ" );
remove_from_PQ(lEdge_to_remove,lData);
}
--mCurrentEdgeCount;
}
--mCurrentEdgeCount;
//the placement is trivial, it's always the point itself
Placement_type lPlacement = lProfile.p0();
vertex_descriptor rResult
= halfedge_collapse_bk_compatibility(lProfile.v0_v1(), Edge_is_constrained_map);
put(Vertex_point_map,rResult,*lPlacement);
Visitor.OnCollapsed(lProfile,rResult);
}
CGAL_SMS_TRACE(0,"Initial edge count: " << mInitialEdgeCount ) ;
}
template<class M,class SP, class VIM, class VPM,class EIM,class ECTM, class CF,class PF,class V>
void EdgeCollapse<M,SP,VIM,VPM,EIM,ECTM,CF,PF,V>::Loop()
{
CGAL_SMS_TRACE(0,"Collapsing edges...") ;
CGAL_SURF_SIMPL_TEST_assertion_code ( size_type lLoop_watchdog = 0 ) ;
CGAL_SURF_SIMPL_TEST_assertion_code ( size_type lNonCollapsableCount = 0 ) ;
//
// Pops and processes each edge from the PQ
//
optional<halfedge_descriptor> lEdge ;
#ifdef CGAL_SURF_SIMPL_INTERMEDIATE_STEPS_PRINTING
int i_rm=0;
#endif
while ( (lEdge = pop_from_PQ()) )
{
CGAL_SURF_SIMPL_TEST_assertion ( lLoop_watchdog ++ < mInitialEdgeCount ) ;
CGAL_SMS_TRACE(1,"Popped " << edge_to_string(*lEdge) ) ;
CGAL_assertion( !is_constrained(*lEdge) );
Profile const& lProfile = create_profile(*lEdge);
Cost_type lCost = get_data(*lEdge).cost();
Visitor.OnSelected(lProfile,lCost,mInitialEdgeCount,mCurrentEdgeCount);
if ( lCost )
{
if ( Should_stop(*lCost,lProfile,mInitialEdgeCount,mCurrentEdgeCount) )
{
Visitor.OnStopConditionReached(lProfile);
CGAL_SMS_TRACE(0,"Stop condition reached with InitialEdgeCount=" << mInitialEdgeCount
<< " CurrentEdgeCount=" << mCurrentEdgeCount
<< " Current Edge: " << edge_to_string(*lEdge)
);
break ;
}
if ( Is_collapse_topologically_valid(lProfile) )
{
// The external function Get_new_vertex_point() is allowed to return an absent point if there is no way to place the vertex
// satisfying its constraints. In that case the remaining vertex is simply left unmoved.
Placement_type lPlacement = get_placement(lProfile);
if ( Is_collapse_geometrically_valid(lProfile,lPlacement) )
{
#ifdef CGAL_SURF_SIMPL_INTERMEDIATE_STEPS_PRINTING
std::cout << "step " << i_rm << " " << source(*lEdge,mSurface)->point() << " " << target(*lEdge,mSurface)->point() << "\n";
#endif
Collapse(lProfile,lPlacement);
#ifdef CGAL_SURF_SIMPL_INTERMEDIATE_STEPS_PRINTING
std::stringstream sstr;
sstr << "debug/P-";
if (i_rm<10) sstr << "0";
if (i_rm<100) sstr << "0";
sstr << i_rm << ".off";
std::ofstream out(sstr.str().c_str());
out << mSurface;
++i_rm;
#endif
}
}
else
{
CGAL_SURF_SIMPL_TEST_assertion_code ( lNonCollapsableCount++ ) ;
Visitor.OnNonCollapsable(lProfile);
CGAL_SMS_TRACE(1,edge_to_string(*lEdge) << " NOT Collapsable" );
}
}
else
{
CGAL_SMS_TRACE(1,edge_to_string(*lEdge) << " uncomputable cost." );
}
}
}
template<class M, class SP, class VIM, class VPM,class EIM,class ECTM, class CF,class PF,class V>
bool EdgeCollapse<M,SP,VIM,VPM,EIM,ECTM,CF,PF,V>::is_border( vertex_descriptor const& aV ) const
{
bool rR = false ;
in_edge_iterator eb, ee ;
for ( boost::tie(eb,ee) = halfedges_around_target(aV,mSurface) ; eb != ee ; ++ eb )
{
halfedge_descriptor lEdge = *eb ;
if ( is_edge_a_border(lEdge) )
{
rR = true ;
break ;
}
}
return rR ;
}
template<class M, class SP, class VIM, class VPM,class EIM,class ECTM, class CF,class PF,class V>
bool EdgeCollapse<M,SP,VIM,VPM,EIM,ECTM,CF,PF,V>::is_border_or_constrained( vertex_descriptor const& aV ) const
{
in_edge_iterator eb, ee ;
for ( boost::tie(eb,ee) = halfedges_around_target(aV,mSurface) ; eb != ee ; ++ eb )
{
halfedge_descriptor lEdge = *eb ;
if ( is_edge_a_border(lEdge) || is_constrained(lEdge) )
return true;
}
return false;
}
template<class M, class SP, class VIM, class VPM,class EIM,class ECTM, class CF,class PF,class V>
bool EdgeCollapse<M,SP,VIM,VPM,EIM,ECTM,CF,PF,V>::is_constrained( vertex_descriptor const& aV ) const
{
in_edge_iterator eb, ee ;
for ( boost::tie(eb,ee) = halfedges_around_target(aV,mSurface) ; eb != ee ; ++ eb )
if ( is_constrained(*eb) ) return true;
return false;
}
// Some edges are NOT collapsable: doing so would break the topological consistency of the mesh.
// This function returns true if a edge 'p->q' can be collapsed.
//
// An edge p->q can be collapsed iff it satisfies the "link condition"
// (as described in the "Mesh Optimization" article of Hoppe et al (1993))
//
// The link condition is as follows: for every vertex 'k' adjacent to both 'p and 'q',
// "p,k,q" is a facet of the mesh.
//
template<class M, class SP, class VIM, class VPM,class EIM,class ECTM, class CF,class PF,class V>
bool EdgeCollapse<M,SP,VIM,VPM,EIM,ECTM,CF,PF,V>::Is_collapse_topologically_valid( Profile const& aProfile )
{
bool rR = true ;
CGAL_SMS_TRACE(3,"Testing topological collapsabilty of p_q=V" << get(Vertex_index_map,aProfile.v0()) << "(%" << degree(aProfile.v0(),mSurface) << ")"
<< "->V" << get(Vertex_index_map,aProfile.v1()) << "(%" << degree(aProfile.v1(),mSurface) << ")"
);
CGAL_SMS_TRACE(4, "is p_q border:" << aProfile.is_v0_v1_a_border() );
CGAL_SMS_TRACE(4, "is q_q border:" << aProfile.is_v1_v0_a_border() );
out_edge_iterator eb1, ee1 ;
out_edge_iterator eb2, ee2 ;
CGAL_SMS_TRACE(4," t=V"
<< ( aProfile.left_face_exists() ? get(Vertex_index_map,aProfile.vL()) : -1 )
<< "(%"
<< ( aProfile.left_face_exists() ? degree(aProfile.vL(),mSurface) : 0 )
<< ")"
);
CGAL_SMS_TRACE(4," b=V"
<< ( aProfile.right_face_exists() ? get(Vertex_index_map,aProfile.vR()) : -1 )
<< "(%"
<< ( aProfile.right_face_exists() ? degree(aProfile.vR(),mSurface) :0 )
<< ")"
);
// Simple tests handling the case of non-manifold situations at a vertex or edge (pinching)
// (even if we advertise one should not use a surface mesh with such features)
if ( aProfile.left_face_exists () )
{
if ( CGAL::is_border( opposite(aProfile.v1_vL(), mSurface), mSurface ) &&
CGAL::is_border( opposite(aProfile.vL_v0(), mSurface), mSurface )
) return false;
if ( aProfile.right_face_exists () &&
CGAL::is_border( opposite(aProfile.vR_v1(), mSurface), mSurface ) &&
CGAL::is_border( opposite(aProfile.v0_vR(), mSurface), mSurface )
) return false;
}
else{
if ( aProfile.right_face_exists () )
{
if ( CGAL::is_border( opposite(aProfile.vR_v1(), mSurface), mSurface ) &&
CGAL::is_border( opposite(aProfile.v0_vR(), mSurface), mSurface )
) return false;
}
else
return false;
}
// The following loop checks the link condition for v0_v1.
// Specifically, that for every vertex 'k' adjacent to both 'p and 'q', 'pkq' is a face of the mesh.
//
for ( boost::tie(eb1,ee1) = halfedges_around_source(aProfile.v0(),mSurface) ; rR && eb1 != ee1 ; ++ eb1 )
{
halfedge_descriptor v0_k = *eb1 ;
if ( v0_k != aProfile.v0_v1() )
{
vertex_descriptor k = target(v0_k,mSurface);
for ( boost::tie(eb2,ee2) = halfedges_around_source(k,mSurface) ; rR && eb2 != ee2 ; ++ eb2 )
{
halfedge_descriptor k_v1 = *eb2 ;
if ( target(k_v1,mSurface) == aProfile.v1() )
{
// At this point we know p-q-k are connected and we need to determine if this triangle is a face of the mesh.
//
// Since the mesh is known to be triangular there are at most two faces sharing the edge p-q.
//
// If p->q is NOT a border edge, the top face is p->q->t where t is target(next(p->q))
// If q->p is NOT a border edge, the bottom face is q->p->b where b is target(next(q->p))
//
// If k is either t or b then p-q-k *might* be a face of the mesh. It won't be if k==t but p->q is border
// or k==b but q->b is a border (because in that case even though there exists triangles p->q->t (or q->p->b)
// they are holes, not faces)
//
bool lIsFace = ( aProfile.vL() == k && aProfile.left_face_exists () )
|| ( aProfile.vR() == k && aProfile.right_face_exists() ) ;
CGAL_SURF_SIMPL_TEST_assertion_code
(
if ( lIsFace )
{
// Is k_v1 the halfedge bounding the face 'k-v1-v0'?
if ( ! is_border(k_v1) && target(next(k_v1,mSurface),mSurface) == aProfile.v0() )
{
CGAL_SURF_SIMPL_TEST_assertion( ! is_border(k_v1) ) ;
CGAL_SURF_SIMPL_TEST_assertion( target(k_v1,mSurface) == aProfile.v1() ) ;
CGAL_SURF_SIMPL_TEST_assertion( target(next(k_v1,mSurface),mSurface) == aProfile.v0() ) ;
CGAL_SURF_SIMPL_TEST_assertion( target(next(next(k_v1,mSurface),mSurface),mSurface) == k ) ;
}
else // or is it the opposite?
{
halfedge_descriptor v1_k = opposite(k_v1,mSurface);
CGAL_SURF_SIMPL_TEST_assertion( ! is_border(v1_k) ) ;
CGAL_SURF_SIMPL_TEST_assertion( target(v1_k,mSurface) == k ) ;
CGAL_SURF_SIMPL_TEST_assertion( target(next(v1_k, mSurface),mSurface) == aProfile.v0() ) ;
CGAL_SURF_SIMPL_TEST_assertion( target(next(next(v1_k, mSurface),mSurface),mSurface) == aProfile.v1() ) ;
}
}
);
if ( !lIsFace )
{
CGAL_SMS_TRACE(3," k=V" << get(Vertex_index_map,k) << " IS NOT in a face with p-q. NON-COLLAPSABLE edge." ) ;
rR = false ;
break ;
}
else
{
CGAL_SMS_TRACE(4," k=V" << get(Vertex_index_map,k) << " is in a face with p-q") ;
}
}
}
}
}
if ( rR )
{
/// ensure two constrained edges cannot get merged
if ( is_edge_adjacent_to_a_constrained_edge(
aProfile, Edge_is_constrained_map) ) return false ;
if ( aProfile.is_v0_v1_a_border() )
{
if ( Is_open_triangle(aProfile.v0_v1()) )
{
rR = false ;
CGAL_SMS_TRACE(3," p-q belongs to an open triangle. NON-COLLAPSABLE edge." ) ;
}
}
else if ( aProfile.is_v1_v0_a_border() )
{
if ( Is_open_triangle(aProfile.v1_v0()) )
{
rR = false ;
CGAL_SMS_TRACE(3," p-q belongs to an open triangle. NON-COLLAPSABLE edge." ) ;
}
}
else
{
if ( is_border(aProfile.v0()) && is_border(aProfile.v1()) )
{
rR = false ;
CGAL_SMS_TRACE(3," both p and q are boundary vertices but p-q is not. NON-COLLAPSABLE edge." ) ;
}
else
{
bool lTetra = Is_tetrahedron(aProfile.v0_v1());
//CGAL_SURF_SIMPL_TEST_assertion( lTetra == mSurface.is_tetrahedron(aProfile.v0_v1()) ) ;
if ( lTetra )
{
rR = false ;
CGAL_SMS_TRACE(3," p-q belongs to a tetrahedron. NON-COLLAPSABLE edge." ) ;
}
if ( next(aProfile.v0_v1(), mSurface) == opposite(prev(aProfile.v1_v0(), mSurface), mSurface) &&
prev(aProfile.v0_v1(), mSurface) == opposite(next(aProfile.v1_v0(), mSurface), mSurface) )
{
CGAL_SMS_TRACE(3," degenerate volume." ) ;
return false ;
}
}
}
}
return rR ;
}
template<class M, class SP, class VIM, class VPM,class EIM,class ECTM, class CF,class PF,class V>
bool EdgeCollapse<M,SP,VIM,VPM,EIM,ECTM,CF,PF,V>::Is_tetrahedron( halfedge_descriptor const& h1 )
{
return is_tetrahedron(h1,mSurface);
}
template<class M, class SP, class VIM, class VPM,class EIM,class ECTM, class CF,class PF,class V>
bool EdgeCollapse<M,SP,VIM,VPM,EIM,ECTM,CF,PF,V>::Is_open_triangle( halfedge_descriptor const& h1 )
{
bool rR = false ;
halfedge_descriptor h2 = next(h1,mSurface);
halfedge_descriptor h3 = next(h2,mSurface);
// First check if it is a triangle
if ( next(h3,mSurface) == h1 )
{
// Now check if it is open
CGAL_SMS_TRACE(4," p-q is a border edge... checking E" << get(Edge_index_map,h2) << " and E" << get(Edge_index_map,h3) ) ;
rR = is_border(h2) && is_border(h3);
CGAL_SURF_SIMPL_TEST_assertion( rR == ( is_border(h1)
&& is_border(next(h1,mSurface))
&& is_border(next(next(h1,mSurface),mSurface))
)
) ;
}
return rR ;
}
// Given triangles 'p0,p1,p2' and 'p0,p2,p3', both shared along edge 'v0-v2',
// determine if they are geometrically valid: that is, the ratio of their
// respective areas is no greater than a max value and the internal
// dihedral angle formed by their supporting planes is no greater than
// a given threshold
template<class M, class SP, class VIM, class VPM,class EIM,class ECTM, class CF,class PF,class V>
bool EdgeCollapse<M,SP,VIM,VPM,EIM,ECTM,CF,PF,V>::are_shared_triangles_valid( Point const& p0, Point const& p1, Point const& p2, Point const& p3 ) const
{
bool rR = false ;
Vector e01 = Traits().construct_vector_3_object()(p0,p1) ;
Vector e02 = Traits().construct_vector_3_object()(p0,p2) ;
Vector e03 = Traits().construct_vector_3_object()(p0,p3) ;
Vector n012 = Traits().construct_cross_product_vector_3_object()(e01,e02);
Vector n023 = Traits().construct_cross_product_vector_3_object()(e02,e03);
FT l012 = Traits().compute_scalar_product_3_object()(n012, n012) ;
FT l023 = Traits().compute_scalar_product_3_object()(n023, n023) ;
FT larger = (std::max)(l012,l023);
FT smaller = (std::min)(l012,l023);
const double cMaxAreaRatio = 1e8 ;
CGAL_SMS_TRACE(4," Testing validity of shared triangles:"
<< "\n p0=" << xyz_to_string(p0) << "\n p1=" << xyz_to_string(p1) << "\n p2=" << xyz_to_string(p2) << "\n p3=" << xyz_to_string(p3)
<< "\n e01=" << xyz_to_string(e01) << "\n e02=" << xyz_to_string(e02) << "\n e03=" << xyz_to_string(e03)
<< "\n n012=" << xyz_to_string(n012) << "\n n023=" << xyz_to_string(n023)
<< "\n l012=" << n_to_string(l012) << "\n l023=" << n_to_string(l023)
);
if ( larger < cMaxAreaRatio * smaller )
{
FT l0123 = Traits().compute_scalar_product_3_object()(n012, n023) ;
CGAL_SMS_TRACE(4,"\n l0123=" << n_to_string(l0123) );
if ( CGAL_NTS is_positive(l0123) )
{
rR = true ;
}
else
{
CGAL_SMS_TRACE(4,"\n lhs: " << n_to_string(( l0123 * l0123 ) / ( l012 * l023 )) << " <= rhs: " << mcMaxDihedralAngleCos2 ) ;
if ( ( l0123 * l0123 ) <= mcMaxDihedralAngleCos2 * ( l012 * l023 ) )
{
rR = true ;
}
}
}
return rR ;
}
// Returns the directed halfedge connecting v0 to v1, if exists.
template<class M, class SP, class VIM, class VPM,class EIM,class ECTM, class CF,class PF,class V>
typename EdgeCollapse<M,SP,VIM,VPM,EIM,ECTM,CF,PF,V>::halfedge_descriptor
EdgeCollapse<M,SP,VIM,VPM,EIM,ECTM,CF,PF,V>::find_connection ( vertex_descriptor const& v0, vertex_descriptor const& v1 ) const
{
out_edge_iterator eb, ee ;
for ( boost::tie(eb,ee) = halfedges_around_source(v0,mSurface) ; eb != ee ; ++ eb )
{
halfedge_descriptor out = *eb ;
if ( target(out,mSurface) == v1 )
return out ;
}
return halfedge_descriptor() ;
}
// Given the edge 'e' around the link for the collapsinge edge "v0-v1", finds the vertex that makes a triangle adjacent to 'e' but exterior to the link (i.e not containing v0 nor v1)
// If 'e' is a null handle OR 'e' is a border edge, there is no such triangle and a null handle is returned.
template<class M, class SP, class VIM, class VPM,class EIM,class ECTM, class CF,class PF,class V>
typename EdgeCollapse<M,SP,VIM,VPM,EIM,ECTM,CF,PF,V>::vertex_descriptor
EdgeCollapse<M,SP,VIM,VPM,EIM,ECTM,CF,PF,V>::find_exterior_link_triangle_3rd_vertex ( halfedge_descriptor const& e, vertex_descriptor const& v0, vertex_descriptor const& v1 ) const
{
vertex_descriptor r ;
if ( handle_assigned(e) )
{
vertex_descriptor ra = target(next(e, mSurface), mSurface);
vertex_descriptor rb = source(prev(e, mSurface), mSurface);
if ( ra == rb && ra != v0 && ra != v1 )
{
r = ra ;
}
else
{
ra = target(next(opposite(e,mSurface), mSurface), mSurface);
rb = source(prev(opposite(e,mSurface), mSurface), mSurface);
if ( ra == rb && ra != v0 && ra != v1 )
{
r = ra ;
}
}
}
return r ;
}
// A collapse is geometrically valid if, in the resulting local mesh no two adjacent triangles form an internal dihedral angle
// greater than a fixed threshold (i.e. triangles do not "fold" into each other)
//
template<class M, class SP, class VIM, class VPM,class EIM,class ECTM, class CF,class PF,class V>
bool EdgeCollapse<M,SP,VIM,VPM,EIM,ECTM,CF,PF,V>::Is_collapse_geometrically_valid( Profile const& aProfile, Placement_type k0)
{
bool rR = true ;
CGAL_SMS_TRACE(3,"Testing geometrical collapsabilty of v0-v1=E" << get(Edge_index_map,aProfile.v0_v1()) );
if ( k0 )
{
// Use the current link to extract all local triangles incident to 'vx' in the collapsed mesh (which at this point doesn't exist yet)
//
typedef typename Profile::vertex_descriptor_vector::const_iterator link_iterator ;
link_iterator linkb = aProfile.link().begin();
link_iterator linke = aProfile.link().end ();
link_iterator linkl = std::prev(linke) ;
for ( link_iterator l = linkb ; l != linke && rR ; ++ l )
{
link_iterator pv = ( l == linkb ? linkl : std::prev (l) );
link_iterator nx = ( l == linkl ? linkb : std::next (l) ) ;
// k0,k1 and k3 are three consecutive vertices along the link.
vertex_descriptor k1 = *pv ;
vertex_descriptor k2 = * l ;
vertex_descriptor k3 = *nx ;
CGAL_SMS_TRACE(4," Screening link vertices k1=V" << get(Vertex_index_map,k1) << " k2=V" << get(Vertex_index_map,k2) << " k3=V" << get(Vertex_index_map,k3) ) ;
halfedge_descriptor e12 = find_connection(k1,k2);
halfedge_descriptor e23 = k3 != k1 ? find_connection(k2,k3) : halfedge_descriptor() ;
// If 'k1-k2-k3' are connected there will be two adjacent triangles 'k0,k1,k2' and 'k0,k2,k3' after the collapse.
if ( handle_assigned(e12) && handle_assigned(e23) )
{
CGAL_SMS_TRACE(4," Link triangles shared" ) ;
if ( !are_shared_triangles_valid( *k0, get_point(k1), get_point(k2), get_point(k3) ) )
{
CGAL_SMS_TRACE(3," Triangles VX-V" << get(Vertex_index_map,k1) << "-V" << get(Vertex_index_map,k2) << " and VX-V" << get(Vertex_index_map,k3) << " are not geometrically valid. Collapse rejected");
rR = false ;
}
}
if ( rR )
{
// Also check the triangles 'k0,k1,k2' and it's adjacent along e12: 'k4,k2,k1', if exist
vertex_descriptor k4 = find_exterior_link_triangle_3rd_vertex(e12,aProfile.v0(),aProfile.v1());
// There is indeed a triangle shared along e12
if ( handle_assigned(k4) )
{
CGAL_SMS_TRACE(4," Found exterior link triangle shared along E" << get(Edge_index_map,e12) << " with third vertex: V" << get(Vertex_index_map,k4) ) ;
if ( !are_shared_triangles_valid( get_point(k1), get_point(k4), get_point(k2), *k0 ) )
{
CGAL_SMS_TRACE(3," Triangles V" << get(Vertex_index_map,k1) << "-V" << get(Vertex_index_map,k4) << " and V" << get(Vertex_index_map,k2) << "-VX are not geometrically valid. Collapse rejected");
rR = false ;
}
}
}
if ( rR )
{
// And finally, check the triangles 'k0,k2,k3' and it's adjacent e23: 'k5,k3,k2' if exist
vertex_descriptor k5 = find_exterior_link_triangle_3rd_vertex(e23,aProfile.v0(),aProfile.v1());
// There is indeed a triangle shared along e12
if ( handle_assigned(k5) )
{
CGAL_SMS_TRACE(4," Found exterior link triangle shared along E" << get(Edge_index_map,e23) << " with third vertex: V" << get(Vertex_index_map,k5) ) ;
if ( !are_shared_triangles_valid( get_point(k2), get_point(k5), get_point(k3), *k0 ) )
{
CGAL_SMS_TRACE(3," Triangles V" << get(Vertex_index_map,k2) << "-V" << get(Vertex_index_map,k5) << " and V" << get(Vertex_index_map,k3) << "-VX are not geometrically valid. Collapse rejected");
rR = false ;
}
}
}
}
}
return rR ;
}
template<class M, class SP, class VIM, class VPM,class EIM,class ECTM, class CF,class PF,class V>
void EdgeCollapse<M,SP,VIM,VPM,EIM,ECTM,CF,PF,V>::Collapse( Profile const& aProfile, Placement_type aPlacement )
{
CGAL_SMS_TRACE(1,"S" << mStep << ". Collapsing " << edge_to_string(aProfile.v0_v1()) ) ;
vertex_descriptor rResult ;
Visitor.OnCollapsing(aProfile,aPlacement);
-- mCurrentEdgeCount ;
CGAL_SURF_SIMPL_TEST_assertion_code( size_type lResultingVertexCount = mSurface.size_of_vertices() ;
size_type lResultingEdgeCount = mSurface.size_of_halfedges() / 2 ;
) ;
// If the top/bottom facets exists, they are removed and the edges v0vt and Q-B along with them.
// In that case their corresponding pairs must be pop off the queue
if ( aProfile.left_face_exists() )
{
halfedge_descriptor lV0VL = primary_edge(aProfile.vL_v0());
if ( is_constrained(lV0VL) ) //make sure a constrained edge will not disappear
lV0VL=primary_edge(aProfile.v1_vL());
CGAL_SMS_TRACE(3,"V0VL E" << get(Edge_index_map,lV0VL)
<< "(V" << get(Vertex_index_map, source(lV0VL,mSurface)) << "->V" << get(Vertex_index_map,target(lV0VL,mSurface)) << ")"
) ;
Edge_data& lData = get_data(lV0VL) ;
if ( lData.is_in_PQ() )
{
CGAL_SMS_TRACE(2,"Removing E" << get(Edge_index_map,lV0VL) << " from PQ" ) ;
remove_from_PQ(lV0VL,lData) ;
}
-- mCurrentEdgeCount ;
CGAL_SURF_SIMPL_TEST_assertion_code( -- lResultingEdgeCount ) ;
}
if ( aProfile.right_face_exists() )
{
halfedge_descriptor lVRV1 = primary_edge(aProfile.vR_v1());
if ( is_constrained(lVRV1) ) //make sure a constrained edge will not disappear
lVRV1=primary_edge(aProfile.v0_vR());
CGAL_SMS_TRACE(3,"V1VRE" << get(Edge_index_map,lVRV1)
<< "(V" << get(Vertex_index_map, source(lVRV1,mSurface)) << "->V" << get(Vertex_index_map, target(lVRV1,mSurface)) << ")"
) ;
Edge_data& lData = get_data(lVRV1) ;
if ( lData.is_in_PQ() )
{
CGAL_SMS_TRACE(2,"Removing E" << get(Edge_index_map,lVRV1) << " from PQ") ;
remove_from_PQ(lVRV1,lData) ;
}
-- mCurrentEdgeCount ;
CGAL_SURF_SIMPL_TEST_assertion_code( -- lResultingEdgeCount ) ;
}
CGAL_SMS_TRACE(1,"Removing:\n v0v1: E" << get(Edge_index_map,aProfile.v0_v1()) << "(V" << get(Vertex_index_map,aProfile.v0()) << "->V" << get(Vertex_index_map,aProfile.v1()) << ")" );
// Perform the actuall collapse.
// This is an external function.
// It's REQUIRED to remove ONLY 1 vertex (P or Q) and edges PQ,PT and QB
// (PT and QB are removed if they are not null).
// All other edges must be kept.
// All directed edges incident to vertex removed are relink to the vertex kept.
rResult = halfedge_collapse_bk_compatibility(aProfile.v0_v1(), Edge_is_constrained_map);
CGAL_SURF_SIMPL_TEST_assertion_code( -- lResultingEdgeCount ) ;
CGAL_SURF_SIMPL_TEST_assertion_code( -- lResultingVertexCount ) ;
CGAL_SURF_SIMPL_TEST_assertion( lResultingEdgeCount * 2 == mSurface.size_of_halfedges() ) ;
CGAL_SURF_SIMPL_TEST_assertion( lResultingVertexCount == mSurface.size_of_vertices() ) ;
CGAL_SURF_SIMPL_TEST_assertion( mSurface.is_valid() && mSurface.is_pure_triangle() ) ;
CGAL_SMS_TRACE(1,"V" << get(Vertex_index_map,rResult) << " kept." ) ;
#ifdef CGAL_SURFACE_SIMPLIFICATION_ENABLE_TRACE
out_edge_iterator eb1, ee1 ;
for ( boost::tie(eb1,ee1) = halfedges_around_source(rResult,mSurface) ; eb1 != ee1 ; ++ eb1 )
CGAL_SMS_TRACE(2, edge_to_string(*eb1) ) ;
#endif
if ( aPlacement )
{
CGAL_SMS_TRACE(1,"New vertex point: " << xyz_to_string(*aPlacement) ) ;
put(Vertex_point_map,rResult,*aPlacement) ;
}
Visitor.OnCollapsed(aProfile,rResult);
Update_neighbors(rResult) ;
CGAL_SMS_DEBUG_CODE ( ++mStep ; )
}
template<class M, class SP, class VIM, class VPM,class EIM,class ECTM, class CF,class PF,class V>
void EdgeCollapse<M,SP,VIM,VPM,EIM,ECTM,CF,PF,V>::Update_neighbors( vertex_descriptor const& aKeptV )
{
CGAL_SMS_TRACE(3,"Updating cost of neighboring edges..." ) ;
//
// (A) Collect all edges to update their cost: all those around each vertex adjacent to the vertex kept
//
typedef std::set<halfedge_descriptor,Compare_id> edges ;
edges lToUpdate(Compare_id(this)) ;
edges lToInsert(Compare_id(this)) ;
// (A.1) Loop around all vertices adjacent to the vertex kept
in_edge_iterator eb1, ee1 ;
for ( boost::tie(eb1,ee1) = halfedges_around_target(aKeptV,mSurface) ; eb1 != ee1 ; ++ eb1 )
{
halfedge_descriptor lEdge1 = *eb1 ;
vertex_descriptor lAdj_k = source(lEdge1,mSurface);
// (A.2) Loop around all edges incident on each adjacent vertex
in_edge_iterator eb2, ee2 ;
for ( boost::tie(eb2,ee2) = halfedges_around_target(lAdj_k,mSurface) ; eb2 != ee2 ; ++ eb2 )
{
halfedge_descriptor lEdge2 = primary_edge(*eb2) ;
Edge_data& lData2 = get_data(lEdge2);
CGAL_SMS_TRACE(4,"Inedge around V" << get(Vertex_index_map,lAdj_k) << edge_to_string(lEdge2) ) ;
// Only edges still in the PQ needs to be updated, the other needs to be re-inserted
if ( lData2.is_in_PQ() )
lToUpdate.insert(lEdge2) ;
else
lToInsert.insert(lEdge2) ;
}
}
//
// (B) Proceed to update the costs.
//
for ( typename edges::iterator it = lToUpdate.begin(), eit = lToUpdate.end() ; it != eit ; ++ it )
{
halfedge_descriptor lEdge = *it;
Edge_data& lData = get_data(lEdge);
Profile const& lProfile = create_profile(lEdge);
lData.cost() = get_cost(lProfile) ;
CGAL_SMS_TRACE(3, edge_to_string(lEdge) << " updated in the PQ") ;
update_in_PQ(lEdge,lData);
}
//
// (C) Insert ignored edges
//
// I think that this should be done for edges eliminated because of the geometric criteria
// and not the topological one.However maintaining such a set might be more expensive
// and hard to be safe ...
for ( typename edges::iterator it = lToInsert.begin(),
eit = lToInsert.end() ; it != eit ; ++ it )
{
halfedge_descriptor lEdge = *it;
if ( is_constrained(lEdge) ) continue; //do not insert constrained edges
Edge_data& lData = get_data(lEdge);
Profile const& lProfile = create_profile(lEdge);
lData.cost() = get_cost(lProfile) ;
CGAL_SMS_TRACE(3, edge_to_string(lEdge) << " re-inserted in the PQ") ;
insert_in_PQ(lEdge,lData);
}
}
} // namespace Surface_mesh_simplification
} //namespace CGAL
#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_DETAIL_EDGE_COLLAPSE_IMPL_H //
// EOF //

53
Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Edge_collapse_visitor_base.h
View File

@ -22,52 +22,37 @@
#include <CGAL/license/Surface_mesh_simplification.h>
#include <CGAL/Surface_mesh_simplification/Detail/Common.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_profile.h>
namespace CGAL {
namespace Surface_mesh_simplification
{
namespace Surface_mesh_simplification {
template<class TM_>
struct Edge_collapse_visitor_base
{
typedef TM_ TM ;
typedef TM_ TM;
typedef Edge_profile<TM> Profile;
typedef boost::graph_traits<TM> GraphTraits;
typedef Edge_profile<TM> Profile ;
typedef boost::graph_traits <TM> GraphTraits ;
typedef typename GraphTraits::edges_size_type size_type ;
typedef typename GraphTraits::vertex_descriptor vertex_descriptor ;
typedef typename GraphTraits::edges_size_type size_type;
typedef typename GraphTraits::vertex_descriptor vertex_descriptor;
typedef typename boost::property_map<TM, CGAL::vertex_point_t>::type Vertex_point_pmap;
typedef typename boost::property_traits<Vertex_point_pmap>::value_type Point;
typedef typename Kernel_traits<Point>::Kernel Kernel ;
typedef typename Kernel::FT FT ;
typedef typename Kernel_traits<Point>::Kernel Kernel;
typedef typename Kernel::FT FT;
void OnStarted( TM& ) {}
void OnFinished ( TM& ) {}
void OnStopConditionReached( Profile const& ) {}
void OnCollected( Profile const&, boost::optional<FT> const& ) {}
void OnSelected( Profile const&, boost::optional<FT> const&, size_type, size_type ) {}
void OnCollapsing(Profile const&, boost::optional<Point> const& ) {}
void OnCollapsed( Profile const&, vertex_descriptor const& ) {}
void OnNonCollapsable(Profile const& ) {}
} ;
void OnStarted(TM&) {}
void OnFinished(TM&) {}
void OnStopConditionReached(const Profile&) {}
void OnCollected(const Profile&, const boost::optional<FT>&) {}
void OnSelected(const Profile&, const boost::optional<FT>&, size_type, size_type) {}
void OnCollapsing(const Profile&, const boost::optional<Point>&) {}
void OnCollapsed(const Profile&, const vertex_descriptor&) {}
void OnNonCollapsable(const Profile&) {}
};
} // namespace Surface_mesh_simplification
} // namespace CGAL
} //namespace CGAL
#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_EDGE_COLLAPSE_VISITOR_BASE_H //
// EOF //
#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_EDGE_COLLAPSE_VISITOR_BASE_H

4
Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/HalfedgeGraph_Polyhedron_3.h
View File

@ -22,11 +22,7 @@
#include <CGAL/license/Surface_mesh_simplification.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/boost/graph/graph_traits_Polyhedron_3.h>
#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_HALFEDGEGRAPH_POLYHEDRON_3_H
// EOF //

46
Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Bounded_normal_change_placement.h
View File

@ -24,56 +24,54 @@
#include <boost/optional.hpp>
namespace CGAL {
namespace Surface_mesh_simplification
{
namespace Surface_mesh_simplification {
template<class Placement>
class Bounded_normal_change_placement
{
public:
typedef typename Placement::TM TM;
typedef typename Placement::TM TM ;
public:
Bounded_normal_change_placement(const Placement& placement = Placement() )
Bounded_normal_change_placement(const Placement& placement = Placement())
: mPlacement(placement)
{}
template <typename Profile>
boost::optional<typename Profile::Point>
operator()( Profile const& aProfile) const
operator()(const Profile& aProfile) const
{
boost::optional<typename Profile::Point> op = mPlacement(aProfile);
if(op){
if(op)
{
// triangles returns the triangles of the star of the vertices of the edge to collapse
// First the two trianges incident to the edge, then the other triangles
// The second vertex of each triangle is the vertex that gets placed
const typename Profile::Triangle_vector& triangles = aProfile.triangles();
if(triangles.size()>2){
if(triangles.size() > 2)
{
typedef typename Profile::Point Point;
typedef typename Profile::Kernel Traits;
typedef typename Traits::Vector_3 Vector;
typename Profile::VertexPointMap ppmap = aProfile.vertex_point_map();
typename Profile::Triangle_vector::const_iterator it = triangles.begin();
if(aProfile.left_face_exists()){
if(aProfile.left_face_exists())
++it;
}
if(aProfile.right_face_exists()){
if(aProfile.right_face_exists())
++it;
}
while(it!= triangles.end()){
while(it!= triangles.end())
{
const typename Profile::Triangle& t = *it;
Point p = get(ppmap,t.v0);
Point q = get(ppmap,t.v1);
Point r = get(ppmap,t.v2);
Point q2 = *op;
Vector eqp = Traits().construct_vector_3_object()(q,p) ;
Vector eqr = Traits().construct_vector_3_object()(q,r) ;
Vector eq2p = Traits().construct_vector_3_object()(q2,p) ;
Vector eq2r = Traits().construct_vector_3_object()(q2,r) ;
Vector eqp = Traits().construct_vector_3_object()(q,p);
Vector eqr = Traits().construct_vector_3_object()(q,r);
Vector eq2p = Traits().construct_vector_3_object()(q2,p);
Vector eq2r = Traits().construct_vector_3_object()(q2,r);
Vector n1 = Traits().construct_cross_product_vector_3_object()(eqp,eqr);
Vector n2 = Traits().construct_cross_product_vector_3_object()(eq2p,eq2r);
@ -88,16 +86,10 @@ public:
}
private:
Placement mPlacement ;
Placement mPlacement;
};
} // namespace Surface_mesh_simplification
} //namespace CGAL
#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_EDGE_COLLAPSE_BOUNDED_NORMAL_CHANGE_PLACEMENT_H

37
Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Constrained_placement.h
View File

@ -22,48 +22,44 @@
#include <CGAL/license/Surface_mesh_simplification.h>
#include <CGAL/Surface_mesh_simplification/Detail/Common.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_profile.h>
namespace CGAL {
namespace Surface_mesh_simplification
{
namespace Surface_mesh_simplification {
template<class BasePlacement, class EdgeIsConstrainedMap>
class Constrained_placement : public BasePlacement
class Constrained_placement
: public BasePlacement
{
public:
EdgeIsConstrainedMap Edge_is_constrained_map;
public:
Constrained_placement(
EdgeIsConstrainedMap map=EdgeIsConstrainedMap(),
BasePlacement base=BasePlacement() )
: BasePlacement(base)
, Edge_is_constrained_map(map)
Constrained_placement(EdgeIsConstrainedMap map=EdgeIsConstrainedMap(),
BasePlacement base=BasePlacement())
: BasePlacement(base),
Edge_is_constrained_map(map)
{}
template <typename Profile>
optional<typename Profile::Point> operator()( Profile const& aProfile ) const
optional<typename Profile::Point> operator()(const Profile& aProfile) const
{
typedef typename Profile::TM TM;
typedef typename CGAL::Halfedge_around_target_iterator<TM> in_edge_iterator;
in_edge_iterator eb, ee ;
for ( boost::tie(eb,ee) = halfedges_around_target(aProfile.v0(),aProfile.surface_mesh());
eb != ee ; ++ eb )
in_edge_iterator eb, ee;
for(boost::tie(eb,ee) = halfedges_around_target(aProfile.v0(),aProfile.surface_mesh());
eb != ee; ++ eb)
{
if( get(Edge_is_constrained_map, edge(*eb,aProfile.surface_mesh())) )
if(get(Edge_is_constrained_map, edge(*eb,aProfile.surface_mesh())))
return get(aProfile.vertex_point_map(),
aProfile.v0());
}
for ( boost::tie(eb,ee) = halfedges_around_target(aProfile.v1(),aProfile.surface_mesh());
eb != ee ; ++ eb )
for(boost::tie(eb,ee) = halfedges_around_target(aProfile.v1(),aProfile.surface_mesh());
eb != ee; ++ eb)
{
if( get(Edge_is_constrained_map, edge(*eb,aProfile.surface_mesh())) )
if(get(Edge_is_constrained_map, edge(*eb,aProfile.surface_mesh())))
return get(aProfile.vertex_point_map(),
aProfile.v1());
}
@ -73,7 +69,6 @@ public:
};
} // namespace Surface_mesh_simplification
} //namespace CGAL
#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_CONSTRAINED_PLACEMENT_H

38
Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Count_ratio_stop_predicate.h
View File

@ -22,14 +22,11 @@
#include <CGAL/license/Surface_mesh_simplification.h>
#include <CGAL/Surface_mesh_simplification/Detail/Common.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_profile.h>
namespace CGAL {
namespace Surface_mesh_simplification
{
namespace Surface_mesh_simplification {
//*******************************************************************************************************************
// -= stopping condition predicate =-
@ -39,42 +36,33 @@ namespace Surface_mesh_simplification
//
//*******************************************************************************************************************
//
// Stops when the ratio of initial to current vertex pairs is below some value.
//
template<class TM_>
class Count_ratio_stop_predicate
{
public:
typedef TM_ TM;
typedef Edge_profile<TM> Profile;
typedef TM_ TM ;
typedef typename boost::graph_traits<TM>::edge_descriptor edge_descriptor;
typedef typename boost::graph_traits<TM>::edges_size_type size_type;
typedef Edge_profile<TM> Profile ;
typedef typename boost::graph_traits<TM>::edge_descriptor edge_descriptor ;
typedef typename boost::graph_traits<TM>::edges_size_type size_type ;
Count_ratio_stop_predicate( double aRatio ) : mRatio(aRatio) {}
Count_ratio_stop_predicate(double aRatio) : mRatio(aRatio) {}
template <typename F>
bool operator()( F const& // aCurrentCost
, Profile const& // aEdgeProfile
, size_type aInitialCount
, size_type aCurrentCount
) const
bool operator()(const F&, // aCurrentCost
const Profile&, // aEdgeProfile
size_type aInitialCount,
size_type aCurrentCount) const
{
return ( static_cast<double>(aCurrentCount) / static_cast<double>(aInitialCount) ) < mRatio ;
return (static_cast<double>(aCurrentCount) / static_cast<double>(aInitialCount)) < mRatio;
}
private:
double mRatio ;
double mRatio;
};
} // namespace Surface_mesh_simplification
} //namespace CGAL
#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_EDGE_COLLAPSE_COUNT_RATIO_STOP_PREDICATE_H //
// EOF //
#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_EDGE_COLLAPSE_COUNT_RATIO_STOP_PREDICATE_H

42
Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Count_stop_predicate.h
View File

@ -18,18 +18,15 @@
// Author(s) : Fernando Cacciola <fernando.cacciola@geometryfactory.com>
//
#ifndef CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_EDGE_COLLAPSE_COUNT_STOP_PREDICATE_H
#define CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_EDGE_COLLAPSE_COUNT_STOP_PREDICATE_H 1
#define CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_EDGE_COLLAPSE_COUNT_STOP_PREDICATE_H
#include <CGAL/license/Surface_mesh_simplification.h>
#include <CGAL/Surface_mesh_simplification/Detail/Common.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_profile.h>
namespace CGAL {
namespace Surface_mesh_simplification
{
namespace Surface_mesh_simplification {
//*******************************************************************************************************************
// -= stopping condition predicate =-
@ -39,45 +36,30 @@ namespace Surface_mesh_simplification
//
//*******************************************************************************************************************
//
// Stops when the number of edges left falls below a given number.
//
template<class TM_>
class Count_stop_predicate
{
public:
typedef TM_ TM;
typedef typename boost::graph_traits<TM>::edges_size_type size_type;
typedef TM_ TM ;
// typedef Edge_profile<TM> Profile ;
typedef typename boost::graph_traits<TM>::edges_size_type size_type ;
// typedef typename Kernel::FT FT ;
public :
Count_stop_predicate( std::size_t aThres ) : mThres(aThres) {}
Count_stop_predicate(std::size_t aThres) : mThres(aThres) {}
template <typename F, typename Profile>
bool operator()( F const& // aCurrentCost
, Profile const& // aEdgeProfile
, std::size_t // aInitialCount
, std::size_t aCurrentCount
) const
bool operator()(const F&, // aCurrentCost
const Profile&, // aEdgeProfile
std::size_t, // aInitialCount
std::size_t aCurrentCount) const
{
return aCurrentCount < mThres ;
return aCurrentCount < mThres;
}
private:
std::size_t mThres ;
std::size_t mThres;
};
} // namespace Surface_mesh_simplification
} //namespace CGAL
#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_EDGE_COLLAPSE_COUNT_STOP_PREDICATE_H //
// EOF //
#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_EDGE_COLLAPSE_COUNT_STOP_PREDICATE_H

752
Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Detail/Lindstrom_Turk_core.h
View File

@ -22,180 +22,740 @@
#include <CGAL/license/Surface_mesh_simplification.h>
#include <vector>
#include <CGAL/Cartesian_converter.h>
#include <CGAL/Surface_mesh_simplification/Detail/Common.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_profile.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/LindstromTurk_params.h>
namespace CGAL {
#include <vector>
//
// This should be in
//
// Implementation of the collapsing cost and placement strategy from:
//
// "Fast and Memory Efficient Polygonal Symplification"
// Peter Lindstrom, Greg Turk
//
namespace Surface_mesh_simplification
{
namespace CGAL {
namespace Surface_mesh_simplification {
template<class TM_, class Profile_>
template<class TM_, class Profile_>
class LindstromTurkCore
{
public:
typedef TM_ TM;
typedef Profile_ Profile;
typedef TM_ TM ;
typedef Profile_ Profile ;
typedef boost::graph_traits<TM> GraphTraits;
typedef boost::graph_traits<TM> GraphTraits ;
typedef typename GraphTraits::vertex_descriptor vertex_descriptor;
typedef typename GraphTraits::halfedge_descriptor halfedge_descriptor;
typedef typename GraphTraits::vertex_descriptor vertex_descriptor ;
typedef typename GraphTraits::halfedge_descriptor halfedge_descriptor ;
typedef LindstromTurk_params Params ;
typedef typename Profile::Point Point ;
typedef LindstromTurk_params Params;
typedef typename Profile::Point Point;
typedef typename Profile::VertexPointMap Vertex_point_pmap;
typedef typename boost::property_traits<Vertex_point_pmap>::value_type TM_Point;
typedef typename Kernel_traits<TM_Point>::Kernel TM_Kernel ;
typedef typename Kernel_traits<TM_Point>::Kernel TM_Kernel;
typedef typename Kernel_traits<Point>::Kernel Kernel;
typedef typename Kernel::Vector_3 Vector ;
typedef typename Kernel::FT FT ;
typedef typename Kernel::Vector_3 Vector;
typedef typename Kernel::FT FT;
typedef optional<FT> Optional_FT ;
typedef optional<Point> Optional_point ;
typedef optional<Vector> Optional_vector ;
typedef optional<FT> Optional_FT;
typedef optional<Point> Optional_point;
typedef optional<Vector> Optional_vector;
typedef MatrixC33<Kernel> Matrix ;
typedef MatrixC33<Kernel> Matrix;
typedef typename Profile::Triangle Triangle ;
typedef typename Profile::vertex_descriptor_vector vertex_descriptor_vector ;
typedef typename Profile::Triangle Triangle;
typedef typename Profile::vertex_descriptor_vector vertex_descriptor_vector;
typedef typename Profile::Triangle_vector ::const_iterator const_triangle_iterator ;
typedef typename Profile::halfedge_descriptor_vector::const_iterator const_border_edge_iterator ;
typedef typename Profile::Triangle_vector::const_iterator const_triangle_iterator;
typedef typename Profile::halfedge_descriptor_vector::const_iterator const_border_edge_iterator;
public:
LindstromTurkCore(Params const& aParams, const Profile& aProfile);
LindstromTurkCore( Params const& aParams, Profile const& aProfile ) ;
Optional_point compute_placement() ;
Optional_FT compute_cost( Optional_point const& p ) ;
Optional_point compute_placement();
Optional_FT compute_cost(const Optional_point& p);
private :
struct Triangle_data
{
Triangle_data( Vector const& aNormalV, FT const& aNormalL ) : NormalV(aNormalV), NormalL(aNormalL) {}
Triangle_data(const Vector& aNormalV, const FT& aNormalL)
: NormalV(aNormalV), NormalL(aNormalL)
{}
Vector NormalV;
FT NormalL;
};
Vector NormalV ;
FT NormalL ;
} ;
struct Boundary_data
{
Boundary_data ( Point s_, Point t_, Vector const& v_, Vector const& n_ ) : s(s_), t(t_), v(v_), n(n_) {}
Boundary_data(Point s_, Point t_, const Vector& v_, const Vector& n_)
: s(s_), t(t_), v(v_), n(n_) {}
Point s, t ;
Vector v, n ;
} ;
typedef std::vector<Triangle_data> Triangle_data_vector ;
typedef std::vector<Boundary_data> Boundary_data_vector ;
Point s, t;
Vector v, n;
};
typedef std::vector<Triangle_data> Triangle_data_vector;
typedef std::vector<Boundary_data> Boundary_data_vector;
private :
void extract_triangle_data();
void extract_boundary_data();
void Extract_triangle_data();
void Extract_boundary_data();
void add_boundary_preservation_constraints(const Boundary_data_vector& aBdry);
void add_volume_preservation_constraints(const Triangle_data_vector& aTriangles);
void add_boundary_and_volume_optimization_constraints(const Boundary_data_vector& aBdry, const Triangle_data_vector& aTriangles);
void add_shape_optimization_constraints(const vertex_descriptor_vector& aLink);
void Add_boundary_preservation_constraints( Boundary_data_vector const& aBdry ) ;
void Add_volume_preservation_constraints( Triangle_data_vector const& aTriangles );
void Add_boundary_and_volume_optimization_constraints( Boundary_data_vector const& aBdry, Triangle_data_vector const& aTriangles ) ;
void Add_shape_optimization_constraints( vertex_descriptor_vector const& aLink ) ;
FT compute_boundary_cost(const Vector& v, const Boundary_data_vector& aBdry);
FT compute_volume_cost(const Vector& v, const Triangle_data_vector& aTriangles);
FT compute_shape_cost(Point const& p, const vertex_descriptor_vector& aLink);
FT Compute_boundary_cost( Vector const& v, Boundary_data_vector const& aBdry ) ;
FT Compute_volume_cost ( Vector const& v, Triangle_data_vector const& aTriangles ) ;
FT Compute_shape_cost ( Point const& p, vertex_descriptor_vector const& aLink ) ;
Point get_point ( vertex_descriptor const& v ) const
Point get_point(const vertex_descriptor v) const
{
return convert(get(mProfile.vertex_point_map(),v));
return convert(get(mProfile.vertex_point_map(), v));
}
static Vector Point_cross_product ( Point const& a, Point const& b )
static Vector point_cross_product(const Point& a, const Point& b)
{
return cross_product(a-ORIGIN,b-ORIGIN);
return cross_product(a-ORIGIN, b-ORIGIN);
}
// This is the (uX)(Xu) product described in the Lindstrom-Turk paper
static Matrix LT_product( Vector const& u )
static Matrix LT_product(const Vector& u)
{
FT a00 = ( u.y()*u.y() ) + ( u.z()*u.z() ) ;
FT a00 = (u.y()*u.y()) + (u.z()*u.z());
FT a01 = -u.x()*u.y();
FT a02 = -u.x()*u.z();
FT a10 = a01 ;
FT a11 = ( u.x()*u.x() ) + ( u.z()*u.z() ) ;
FT a10 = a01;
FT a11 = (u.x()*u.x()) + (u.z()*u.z());
FT a12 = - u.y() * u.z();
FT a20 = a02 ;
FT a21 = a12 ;
FT a22 = ( u.x()*u.x() ) + ( u.y()*u.y() ) ;
FT a20 = a02;
FT a21 = a12;
FT a22 = (u.x()*u.x()) + (u.y()*u.y());
return Matrix(a00,a01,a02
,a10,a11,a12
,a20,a21,a22
);
return Matrix(a00, a01, a02,
a10, a11, a12,
a20, a21, a22);
}
static FT big_value() { return static_cast<FT>((std::numeric_limits<double>::max)()) ; }
static FT big_value() { return static_cast<FT>((std::numeric_limits<double>::max)()); }
static bool is_finite ( FT const& n ) { return CGAL_NTS is_finite(n) ; }
static bool is_finite ( Point const& p ) { return is_finite(p.x()) && is_finite(p.y()) && is_finite(p.z()) ; }
static bool is_finite ( Vector const& v ) { return is_finite(v.x()) && is_finite(v.y()) && is_finite(v.z()) ; }
static bool is_finite ( Matrix const& m ) { return is_finite(m.r0()) && is_finite(m.r1()) && is_finite(m.r2()) ; }
static bool is_finite(const FT& n) { return CGAL_NTS is_finite(n); }
static bool is_finite(const Point& p) { return is_finite(p.x()) && is_finite(p.y()) && is_finite(p.z()); }
static bool is_finite(const Vector& v) { return is_finite(v.x()) && is_finite(v.y()) && is_finite(v.z()); }
static bool is_finite(const Matrix& m) { return is_finite(m.r0()) && is_finite(m.r1()) && is_finite(m.r2()); }
template<class T>
static optional<T> filter_infinity ( T const& n ) { return is_finite(n) ? optional<T>(n) : optional<T>() ; }
static optional<T> filter_infinity (T const& n) { return is_finite(n) ? optional<T>(n) : optional<T>(); }
TM& surface() const { return mProfile.surface() ; }
TM& surface() const { return mProfile.surface(); }
private:
Params const& mParams;
const Profile& mProfile;
Params const& mParams ;
Profile const& mProfile ;
void Add_constraint_if_alpha_compatible( Vector const& Ai, FT const& bi ) ;
void Add_constraint_from_gradient ( Matrix const& H, Vector const& c ) ;
void add_constraint_if_alpha_compatible(const Vector& Ai, const FT& bi);
void add_constraint_from_gradient(const Matrix& H, const Vector& c);
private:
Triangle_data_vector mTriangle_data;
Boundary_data_vector mBdry_data;
Triangle_data_vector mTriangle_data ;
Boundary_data_vector mBdry_data ;
int mConstraints_n;
Matrix mConstraints_A;
Vector mConstraints_b;
int mConstraints_n ;
Matrix mConstraints_A ;
Vector mConstraints_b ;
Cartesian_converter<TM_Kernel,Kernel> convert ;
Cartesian_converter<TM_Kernel, Kernel> convert;
FT mSquared_cos_alpha;
FT mSquared_sin_alpha;
};
template<class TM, class K>
LindstromTurkCore<TM,K>::
LindstromTurkCore(Params const& aParams, const Profile& aProfile)
:
mParams(aParams),
mProfile(aProfile),
mConstraints_n(0),
mConstraints_A(NULL_MATRIX),
mConstraints_b(NULL_VECTOR)
{
double alpha = 1.0 * CGAL_PI / 180.0;
FT cos_alpha = std::cos(alpha);
FT sin_alpha = std::sin(alpha);
mSquared_cos_alpha = cos_alpha * cos_alpha;
mSquared_sin_alpha = sin_alpha * sin_alpha;
extract_triangle_data();
extract_boundary_data();
}
template<class TM, class K>
void
LindstromTurkCore<TM,K>::
extract_boundary_data()
{
mBdry_data.reserve(mProfile.border_edges().size());
for(const_border_edge_iterator it = mProfile.border_edges().begin(), eit = mProfile.border_edges().end(); it != eit; ++ it)
{
halfedge_descriptor border_edge = *it;
halfedge_descriptor face_edge = opposite(border_edge,surface());
vertex_descriptor sv = source(face_edge,surface());
vertex_descriptor tv = target(face_edge,surface());
const Point& sp = get_point(sv);
const Point& tp = get_point(tv);
Vector v = tp - sp;
Vector n = point_cross_product(tp,sp);
CGAL_SMS_LT_TRACE(1, " Boundary edge. S:" << xyz_to_string(sp) << " T:" << xyz_to_string(tp)
<< " V:" << xyz_to_string(v) << " N:" << xyz_to_string(n));
mBdry_data.push_back(Boundary_data(sp,tp,v,n));
}
}
template<class TM, class K>
void
LindstromTurkCore<TM,K>::
extract_triangle_data()
{
mTriangle_data.reserve(mProfile.triangles().size());
for(const_triangle_iterator it = mProfile.triangles().begin(), eit = mProfile.triangles().end(); it != eit; ++ it)
{
Triangle const& tri = *it;
const Point& p0 = get_point(tri.v0);
const Point& p1 = get_point(tri.v1);
const Point& p2 = get_point(tri.v2);
Vector v01 = p1 - p0;
Vector v02 = p2 - p0;
Vector lNormalV = cross_product(v01,v02);
FT lNormalL = point_cross_product(p0,p1) * (p2-ORIGIN);
CGAL_SMS_LT_TRACE(1, " Extracting triangle v" << tri.v0 << "->v" << tri.v1 << "->v" << tri.v2
<< " N:" << xyz_to_string(lNormalV) << " L:" << n_to_string(lNormalL));
mTriangle_data.push_back(Triangle_data(lNormalV,lNormalL));
}
}
template<class TM, class K>
typename LindstromTurkCore<TM,K>::Optional_point
LindstromTurkCore<TM,K>::
compute_placement()
{
Optional_point rPlacementP;
Optional_vector lPlacementV;
CGAL_SMS_LT_TRACE(0, "Computing LT placement for E" << mProfile.v0_v1()
<< " (V" << mProfile.v0() << "->V" << mProfile.v1() << ")");
// Each vertex constraint is an equation of the form: Ai * v = bi
// Where 'v' is a vector representing the vertex,
// 'Ai' is a (row) vector and 'bi' a scalar.
//
// The vertex is completely determined with 3 such constraints,
// so is the solution to the folloing system:
//
// A.r0(). * v = b0
// A1 * v = b1
// A2 * v = b2
//
// Which in matrix form is : A * v = b
//
// (with 'A' a 3x3 matrix and 'b' a vector)
//
// The member variable mConstrinas contains A and b. Indidivual constraints (Ai,bi) can be added to it.
// Once 3 such constraints have been added 'v' is directly solved a:
//
// v = b*inverse(A)
//
// A constraint (Ai,bi) must be alpha-compatible with the previously added constraints (see Paper); if it's not, is discarded.
//
if(mBdry_data.size() > 0)
add_boundary_preservation_constraints(mBdry_data);
if(mConstraints_n < 3)
add_volume_preservation_constraints(mTriangle_data);
if(mConstraints_n < 3)
add_boundary_and_volume_optimization_constraints(mBdry_data,mTriangle_data);
if(mConstraints_n < 3)
add_shape_optimization_constraints(mProfile.link());
// It might happen that there were not enough alpha-compatible constraints.
// In that case there is simply no good vertex placement
if(mConstraints_n == 3)
{
// If the matrix is singular it's inverse cannot be computed so an 'absent' value is returned.
optional<Matrix> lOptional_Ai = inverse_matrix(mConstraints_A);
if(lOptional_Ai)
{
const Matrix& lAi = *lOptional_Ai;
CGAL_SMS_LT_TRACE(2," b: " << xyz_to_string(mConstraints_b));
CGAL_SMS_LT_TRACE(2," inv(A): " << matrix_to_string(lAi));
lPlacementV = filter_infinity(mConstraints_b * lAi);
CGAL_SMS_LT_TRACE(0," New vertex point: " << xyz_to_string(*lPlacementV));
}
else
{
CGAL_SMS_LT_TRACE(1," Can't solve optimization, singular system.");
}
}
else
{
CGAL_SMS_LT_TRACE(1," Can't solve optimization, not enough alpha-compatible constraints.");
}
if(lPlacementV)
rPlacementP = ORIGIN + (*lPlacementV);
return rPlacementP;
}
template<class TM, class K>
typename LindstromTurkCore<TM,K>::Optional_FT
LindstromTurkCore<TM,K>::
compute_cost(const Optional_point& aP)
{
Optional_FT rCost;
if(aP)
{
CGAL_SMS_LT_TRACE(0,"Computing LT cost for E" << mProfile.v0_v1());
Vector lV = (*aP) - ORIGIN;
FT lSquaredLength = squared_distance(mProfile.p0(),mProfile.p1());
FT lBdryCost = compute_boundary_cost(lV ,mBdry_data);
FT lVolumeCost = compute_volume_cost(lV ,mTriangle_data);
FT lShapeCost = compute_shape_cost(*aP,mProfile.link());
FT lTotalCost = FT(mParams.VolumeWeight) * lVolumeCost
+ FT(mParams.BoundaryWeight) * lBdryCost * lSquaredLength
+ FT(mParams.ShapeWeight) * lShapeCost * lSquaredLength * lSquaredLength;
rCost = filter_infinity(lTotalCost);
CGAL_SMS_LT_TRACE(0, " Squared edge length: " << n_to_string(lSquaredLength)
<< "\n Boundary cost: " << n_to_string(lBdryCost) << " weight: " << mParams.BoundaryWeight
<< "\n Volume cost: " << n_to_string(lVolumeCost) << " weight: " << mParams.VolumeWeight
<< "\n Shape cost: " << n_to_string(lShapeCost) << " weight: " << mParams.ShapeWeight
<< "\n TOTAL COST: " << n_to_string(lTotalCost));
}
return rCost;
}
template<class TM, class K>
void
LindstromTurkCore<TM,K>::
add_boundary_preservation_constraints(const Boundary_data_vector& aBdry)
{
if(aBdry.size() > 0)
{
Vector e1 = NULL_VECTOR;
Vector e2 = NULL_VECTOR;
FT e1x = FT(0),
e1y = FT(0),
e1z = FT(0);
for(typename Boundary_data_vector::const_iterator it = aBdry.begin(); it != aBdry.end(); ++ it)
{
e1 = e1 + it->v;
e2 = e2 + it->n;
FT vx = it->v.x();
FT vy = it->v.y();
FT vz = it->v.z();
e1x = e1x + vx;
e1y = e1y + vy;
e1z = e1z + vz;
CGAL_SMS_LT_TRACE(1," vx:" << n_to_string(vx) << " vy:" << n_to_string(vy) << " vz:" << n_to_string(vz) << " e1x:"
<< n_to_string(e1x) << " e1y:" << n_to_string(e1y) << " e1z:" << n_to_string(e1z));
}
Matrix H = LT_product(e1);
Vector c = cross_product(e1,e2);
CGAL_SMS_LT_TRACE(1, " Adding boundary preservation constraint."
<< "\n SumV:" << xyz_to_string(e1)
<< "\n SumN:" << xyz_to_string(e2)
<< "\n H:" << matrix_to_string(H)
<< "\n c:" << xyz_to_string(c));
add_constraint_from_gradient(H,c);
}
}
template<class TM, class K>
void
LindstromTurkCore<TM,K>::
add_volume_preservation_constraints(const Triangle_data_vector& aTriangles)
{
CGAL_SMS_LT_TRACE(1," Adding volume preservation constraint. " << aTriangles.size() << " triangles.");
Vector lSumV = NULL_VECTOR;
FT lSumL(0);
for(typename Triangle_data_vector::const_iterator it = aTriangles.begin(), eit = aTriangles.end(); it != eit; ++it)
{
lSumV = lSumV + it->NormalV;
lSumL = lSumL + it->NormalL;
}
CGAL_SMS_LT_TRACE(1, " SumV:" << xyz_to_string(lSumV) << "\n SumL:" << n_to_string(lSumL));
add_constraint_if_alpha_compatible(lSumV,lSumL);
}
template<class TM, class K>
void
LindstromTurkCore<TM,K>::
add_boundary_and_volume_optimization_constraints(const Boundary_data_vector& aBdry,
const Triangle_data_vector& aTriangles)
{
CGAL_SMS_LT_TRACE(1," Adding boundary and volume optimization constraints. ");
Matrix H = NULL_MATRIX;
Vector c = NULL_VECTOR;
// Volume optimization
for(typename Triangle_data_vector::const_iterator it = aTriangles.begin(), eit = aTriangles.end(); it != eit; ++it)
{
Triangle_data const& lTri = *it;
H += direct_product(lTri.NormalV, lTri.NormalV);
c = c - (lTri.NormalL * lTri.NormalV);
}
CGAL_SMS_LT_TRACE(2," Hv:" << matrix_to_string(H) << "\n cv:" << xyz_to_string(c));
if(aBdry.size() > 0)
{
// Boundary optimization
Matrix Hb = NULL_MATRIX;
Vector cb = NULL_VECTOR;
for(typename Boundary_data_vector::const_iterator it = aBdry.begin(); it != aBdry.end(); ++ it)
{
Matrix H = LT_product(it->v);
Vector c = cross_product(it->v,it->n);
Hb += H;
cb = cb + c;
}
CGAL_SMS_LT_TRACE(2," Hb:" << matrix_to_string(Hb) << "\n cb:" << xyz_to_string(cb));
// Weighted average
FT lScaledBoundaryWeight = FT(9) * FT(mParams.BoundaryWeight) * squared_distance(mProfile.p0(),mProfile.p1()) ;
H *= mParams.VolumeWeight;
c = c * mParams.VolumeWeight;
H += lScaledBoundaryWeight * Hb;
c = c + (lScaledBoundaryWeight * cb);
CGAL_SMS_LT_TRACE(2," H:" << matrix_to_string(H) << "\n c:" << xyz_to_string(c));
CGAL_SMS_LT_TRACE(2," VolW:" << mParams.VolumeWeight << " BdryW:" << mParams.BoundaryWeight << " ScaledBdryW:" << lScaledBoundaryWeight);
}
add_constraint_from_gradient(H,c);
}
template<class TM, class K>
void
LindstromTurkCore<TM,K>::
add_shape_optimization_constraints(const vertex_descriptor_vector& aLink)
{
FT s(double(aLink.size()));
Matrix H (s, 0.0, 0.0,
0.0, s, 0.0,
0.0, 0.0, s);
Vector c = NULL_VECTOR;
for(typename vertex_descriptor_vector::const_iterator it = aLink.begin(), eit = aLink.end(); it != eit; ++it)
c = c + (ORIGIN - get_point(*it));
CGAL_SMS_LT_TRACE(1," Adding shape optimization constraint. Shape vector: " << xyz_to_string(c));
add_constraint_from_gradient(H,c);
}
template<class TM, class K>
typename LindstromTurkCore<TM,K>::FT
LindstromTurkCore<TM,K>::
compute_boundary_cost(const Vector& v,
const Boundary_data_vector& aBdry)
{
FT rCost(0);
for(typename Boundary_data_vector::const_iterator it = aBdry.begin(); it != aBdry.end(); ++ it)
{
Vector u = (it->t - ORIGIN) - v;
Vector c = cross_product(it->v,u);
rCost += c*c;
}
return rCost / FT(4);
}
template<class TM, class K>
typename LindstromTurkCore<TM,K>::FT
LindstromTurkCore<TM,K>::
compute_volume_cost(const Vector& v,
const Triangle_data_vector& aTriangles)
{
FT rCost(0);
for(typename Triangle_data_vector::const_iterator it = aTriangles.begin(), eit = aTriangles.end(); it != eit; ++it)
{
Triangle_data const& lTri = *it;
FT lF = lTri.NormalV * v - lTri.NormalL;
rCost += (lF * lF);
}
return rCost / FT(36);
}
template<class TM, class K>
typename LindstromTurkCore<TM,K>::FT
LindstromTurkCore<TM,K>::
compute_shape_cost(const Point& p,
const vertex_descriptor_vector& aLink)
{
FT rCost(0);
for(typename vertex_descriptor_vector::const_iterator it = aLink.begin(), eit = aLink.end(); it != eit; ++it)
rCost += squared_distance(p,get_point(*it));
return rCost;
}
template<class TM, class K>
void
LindstromTurkCore<TM,K>::
add_constraint_if_alpha_compatible(const Vector& Ai,
const FT& bi)
{
CGAL_SMS_LT_TRACE(3, " Adding new constraints if alpha-compatible.\n Ai: " << xyz_to_string(Ai) << "\n bi:" << n_to_string(bi) << ")");
FT slai = Ai*Ai;
CGAL_SMS_LT_TRACE(3, "\n slai: " << n_to_string(slai) << ")");
if(is_finite(slai))
{
FT l = CGAL_NTS sqrt(slai);
CGAL_SMS_LT_TRACE(3, " l: " << n_to_string(l));
if(!CGAL_NTS is_zero(l))
{
Vector Ain = Ai / l;
FT bin = bi / l;
CGAL_SMS_LT_TRACE(3," Ain: " << xyz_to_string(Ain) << " bin:" << n_to_string(bin));
bool lAddIt = true;
if(mConstraints_n == 1)
{
FT d01 = mConstraints_A.r0() * Ai ;
FT sla0 = mConstraints_A.r0() * mConstraints_A.r0();
FT sd01 = d01 * d01;
FT max = sla0 * slai * mSquared_cos_alpha;
CGAL_SMS_LT_TRACE(3," Second constraint. d01: " << n_to_string(d01) << " sla0:" << n_to_string(sla0) << " sd01:" << n_to_string(sd01) << " max:" << n_to_string(max));
if(sd01 > max)
lAddIt = false;
}
else if(mConstraints_n == 2)
{
Vector N = cross_product(mConstraints_A.r0(),mConstraints_A.r1());
FT dc012 = N * Ai;
FT slc01 = N * N;
FT sdc012 = dc012 * dc012;
FT min = slc01 * slai * mSquared_sin_alpha;
CGAL_SMS_LT_TRACE(3, " Third constraint. N: " << xyz_to_string(N) << " dc012:" << n_to_string(dc012) << " slc01:" << n_to_string(slc01)
<< " sdc012:" << n_to_string(sdc012) << " min:" << n_to_string(min));
if(sdc012 <= min)
lAddIt = false;
}
if(lAddIt)
{
switch(mConstraints_n)
{
case 0:
mConstraints_A.r0() = Ain;
mConstraints_b = Vector(bin,mConstraints_b.y(),mConstraints_b.z());
break;
case 1:
mConstraints_A.r1() = Ain;
mConstraints_b = Vector(mConstraints_b.x(),bin,mConstraints_b.z());
break;
case 2:
mConstraints_A.r2() = Ain;
mConstraints_b = Vector(mConstraints_b.x(),mConstraints_b.y(),bin);
break;
}
CGAL_SMS_LT_TRACE(3," Accepting # " << mConstraints_n << " A:" << matrix_to_string(mConstraints_A) << " b:" << xyz_to_string(mConstraints_b));
++ mConstraints_n;
}
else
{
CGAL_SMS_LT_TRACE(3," INCOMPATIBLE. Discarded");
}
}
else
{
CGAL_SMS_LT_TRACE(3," l is ZERO. Discarded");
}
}
else
{
CGAL_SMS_LT_TRACE(3," OVERFLOW. Discarded");
}
}
template<class V>
int index_of_max_component (V const& v)
{
typedef typename Kernel_traits<V>::Kernel::FT FT;
int i = 0;
FT max = v.x();
if(max < v.y())
{
max = v.y();
i = 1;
}
if(max < v.z())
{
max = v.z();
i = 2;
}
return i;
}
template<class TM, class K>
void
LindstromTurkCore<TM,K>::
add_constraint_from_gradient(const Matrix& H,
const Vector& c)
{
CGAL_SMS_LT_TRACE(3," Adding constraint from gradient. Current n=" << mConstraints_n);
CGAL_precondition(mConstraints_n >= 0 && mConstraints_n<=2);
switch(mConstraints_n)
{
case 0:
add_constraint_if_alpha_compatible(H.r0(),-c.x());
add_constraint_if_alpha_compatible(H.r1(),-c.y());
add_constraint_if_alpha_compatible(H.r2(),-c.z());
break;
case 1:
{
const Vector& A0 = mConstraints_A.r0();
CGAL_assertion(A0 != NULL_VECTOR);
Vector AbsA0(CGAL::abs(A0.x()),
CGAL::abs(A0.y()),
CGAL::abs(A0.z()));
Vector Q0;
switch(index_of_max_component(AbsA0))
{
// Since A0 is guaranteed to be non-zero, the denominators here are known to be non-zero too.
case 0:
Q0 = Vector(- A0.z()/A0.x(), 0, 1);
break;
case 1:
Q0 = Vector(0,- A0.z()/A0.y(), 1);
break;
case 2:
Q0 = Vector(1, 0, - A0.x()/A0.z());
break;
default:
Q0 = NULL_VECTOR; // This should never happen!
}
CGAL_SMS_LT_TRACE(3," Q0:" << xyz_to_string(Q0));
CGAL_assertion(Q0 != NULL_VECTOR);
Vector Q1 = cross_product(A0,Q0);
Vector A1 = H * Q0;
Vector A2 = H * Q1;
FT b1 = - (Q0 * c);
FT b2 = - (Q1 * c);
CGAL_SMS_LT_TRACE(3, " Q1:" << xyz_to_string(Q1) << "\n A1: " << xyz_to_string(A1) << "\n A2:" << xyz_to_string(A2)
<< "\n b1:" << n_to_string(b1) << "\n b2:" << n_to_string(b2));
add_constraint_if_alpha_compatible(A1,b1);
add_constraint_if_alpha_compatible(A2,b2);
}
break;
case 2:
{
Vector Q = cross_product(mConstraints_A.r0(),mConstraints_A.r1());
Vector A2 = H * Q;
FT b2 = - (Q * c);
CGAL_SMS_LT_TRACE(3," Q:" << xyz_to_string(Q) << "\n A2: " << xyz_to_string(A2) << "\n b2:" << n_to_string(b2));
add_constraint_if_alpha_compatible(A2,b2);
}
break;
}
}
} // namespace Surface_mesh_simplification
} // namespace CGAL
} //namespace CGAL
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Detail/Lindstrom_Turk_core_impl.h>
#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_EDGE_COLLAPSE_DETAIL_LINDSTROM_TURK_CORE_H //
// EOF //
#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_EDGE_COLLAPSE_DETAIL_LINDSTROM_TURK_CORE_H

620
Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Detail/Lindstrom_Turk_core_impl.h
View File

@ -1,620 +0,0 @@
// Copyright (c) 2005, 2006 Fernando Luis Cacciola Carballal. All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
// You can redistribute it and/or modify it under the terms of the GNU
// General Public License as published by the Free Software Foundation,
// either version 3 of the License, or (at your option) any later version.
//
// Licensees holding a valid Surface_mesh_simplification license may use this file in
// accordance with the Surface_mesh_simplification license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0+
//
// Author(s) : Fernando Cacciola <fernando_cacciola@ciudad.com.ar>
//
#ifndef CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_EDGE_COLLAPSE_DETAIL_LINDSTROM_TURK_CORE_IMPL_H
#define CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_EDGE_COLLAPSE_DETAIL_LINDSTROM_TURK_CORE_IMPL_H 1
#include <CGAL/license/Surface_mesh_simplification.h>
namespace CGAL {
//
// Implementation of the strategy from:
//
// "Fast and Memory Efficient Polygonal Symplification"
// Peter Lindstrom, Greg Turk
//
namespace Surface_mesh_simplification
{
template<class TM, class K>
LindstromTurkCore<TM,K>::LindstromTurkCore( Params const& aParams, Profile const& aProfile )
:
mParams(aParams)
,mProfile(aProfile)
,mConstraints_n(0)
,mConstraints_A(NULL_MATRIX)
,mConstraints_b(NULL_VECTOR)
{
double alpha = 1.0 * CGAL_PI / 180.0;
FT cos_alpha = std::cos(alpha);
FT sin_alpha = std::sin(alpha);
mSquared_cos_alpha = cos_alpha * cos_alpha ;
mSquared_sin_alpha = sin_alpha * sin_alpha ;
Extract_triangle_data();
Extract_boundary_data();
}
template<class TM, class K>
void LindstromTurkCore<TM,K>::Extract_boundary_data()
{
mBdry_data.reserve(mProfile.border_edges().size());
for ( const_border_edge_iterator it = mProfile.border_edges().begin(), eit = mProfile.border_edges().end() ; it != eit ; ++ it )
{
halfedge_descriptor border_edge = *it ;
halfedge_descriptor face_edge = opposite(border_edge,surface()) ;
vertex_descriptor sv = source(face_edge,surface());
vertex_descriptor tv = target(face_edge,surface());
Point const& sp = get_point(sv);
Point const& tp = get_point(tv);
Vector v = tp - sp ;
Vector n = Point_cross_product(tp,sp) ;
CGAL_SMS_LT_TRACE(1," Boundary edge. S:" << xyz_to_string(sp) << " T:" << xyz_to_string(tp)
<< " V:" << xyz_to_string(v) << " N:" << xyz_to_string(n)
) ;
mBdry_data.push_back( Boundary_data(sp,tp,v,n) ) ;
}
}
template<class TM, class K>
void LindstromTurkCore<TM,K>::Extract_triangle_data()
{
mTriangle_data.reserve(mProfile.triangles().size());
for ( const_triangle_iterator it = mProfile.triangles().begin(), eit = mProfile.triangles().end() ; it != eit ; ++ it )
{
Triangle const& tri = *it ;
Point const& p0 = get_point(tri.v0);
Point const& p1 = get_point(tri.v1);
Point const& p2 = get_point(tri.v2);
Vector v01 = p1 - p0 ;
Vector v02 = p2 - p0 ;
Vector lNormalV = cross_product(v01,v02);
FT lNormalL = Point_cross_product(p0,p1) * (p2-ORIGIN);
CGAL_SMS_LT_TRACE(1," Extracting triangle v" << tri.v0 << "->v" << tri.v1 << "->v" << tri.v2
<< " N:" << xyz_to_string(lNormalV) << " L:" << n_to_string(lNormalL)
);
mTriangle_data.push_back(Triangle_data(lNormalV,lNormalL));
}
}
template<class TM, class K>
typename LindstromTurkCore<TM,K>::Optional_point LindstromTurkCore<TM,K>::compute_placement()
{
Optional_point rPlacementP ;
Optional_vector lPlacementV ;
CGAL_SMS_LT_TRACE(0,"Computing LT placement for E" << mProfile.v0_v1() << " (V" << mProfile.v0() << "->V" << mProfile.v1() << ")" );
//
// Each vertex constraint is an equation of the form: Ai * v = bi
// Where 'v' is a vector representing the vertex,
// 'Ai' is a (row) vector and 'bi' a scalar.
//
// The vertex is completely determined with 3 such constraints,
// so is the solution to the folloing system:
//
// A.r0(). * v = b0
// A1 * v = b1
// A2 * v = b2
//
// Which in matrix form is : A * v = b
//
// (with 'A' a 3x3 matrix and 'b' a vector)
//
// The member variable mConstrinas contains A and b. Indidivual constraints (Ai,bi) can be added to it.
// Once 3 such constraints have been added 'v' is directly solved a:
//
// v = b*inverse(A)
//
// A constraint (Ai,bi) must be alpha-compatible with the previously added constraints (see Paper); if it's not, is discarded.
//
if ( mBdry_data.size() > 0 )
Add_boundary_preservation_constraints(mBdry_data);
if ( mConstraints_n < 3 )
Add_volume_preservation_constraints(mTriangle_data);
if ( mConstraints_n < 3 )
Add_boundary_and_volume_optimization_constraints(mBdry_data,mTriangle_data);
if ( mConstraints_n < 3 )
Add_shape_optimization_constraints(mProfile.link());
// It might happen that there were not enough alpha-compatible constraints.
// In that case there is simply no good vertex placement
if ( mConstraints_n == 3 )
{
// If the matrix is singular it's inverse cannot be computed so an 'absent' value is returned.
optional<Matrix> lOptional_Ai = inverse_matrix(mConstraints_A);
if ( lOptional_Ai )
{
Matrix const& lAi = *lOptional_Ai ;
CGAL_SMS_LT_TRACE(2," b: " << xyz_to_string(mConstraints_b) );
CGAL_SMS_LT_TRACE(2," inv(A): " << matrix_to_string(lAi) );
lPlacementV = filter_infinity(mConstraints_b * lAi) ;
CGAL_SMS_LT_TRACE(0," New vertex point: " << xyz_to_string(*lPlacementV) );
}
else
{
CGAL_SMS_LT_TRACE(1," Can't solve optimization, singular system.");
}
}
else
{
CGAL_SMS_LT_TRACE(1," Can't solve optimization, not enough alpha-compatible constraints.");
}
if ( lPlacementV )
rPlacementP = ORIGIN + (*lPlacementV) ;
return rPlacementP;
}
template<class TM, class K>
typename LindstromTurkCore<TM,K>::Optional_FT LindstromTurkCore<TM,K>::compute_cost( Optional_point const& aP )
{
Optional_FT rCost ;
if ( aP )
{
CGAL_SMS_LT_TRACE(0,"Computing LT cost for E" << mProfile.v0_v1() );
Vector lV = (*aP) - ORIGIN ;
FT lSquaredLength = squared_distance(mProfile.p0(),mProfile.p1());
FT lBdryCost = Compute_boundary_cost(lV ,mBdry_data);
FT lVolumeCost = Compute_volume_cost (lV ,mTriangle_data);
FT lShapeCost = Compute_shape_cost (*aP,mProfile.link());
FT lTotalCost = FT(mParams.VolumeWeight) * lVolumeCost
+ FT(mParams.BoundaryWeight) * lBdryCost * lSquaredLength
+ FT(mParams.ShapeWeight) * lShapeCost * lSquaredLength * lSquaredLength ;
rCost = filter_infinity(lTotalCost);
CGAL_SMS_LT_TRACE(0, " Squared edge length: " << n_to_string(lSquaredLength)
<< "\n Boundary cost: " << n_to_string(lBdryCost) << " weight: " << mParams.BoundaryWeight
<< "\n Volume cost: " << n_to_string(lVolumeCost) << " weight: " << mParams.VolumeWeight
<< "\n Shape cost: " << n_to_string(lShapeCost) << " weight: " << mParams.ShapeWeight
<< "\n TOTAL COST: " << n_to_string(lTotalCost)
);
}
return rCost;
}
template<class TM, class K>
void LindstromTurkCore<TM,K>::Add_boundary_preservation_constraints( Boundary_data_vector const& aBdry )
{
if ( aBdry.size() > 0 )
{
Vector e1 = NULL_VECTOR ;
Vector e2 = NULL_VECTOR ;
FT e1x = FT(0.0), e1y = FT(0.0), e1z = FT(0.0) ;
for ( typename Boundary_data_vector::const_iterator it = aBdry.begin() ; it != aBdry.end() ; ++ it )
{
e1 = e1 + it->v ;
e2 = e2 + it->n ;
FT vx = it->v.x() ;
FT vy = it->v.y() ;
FT vz = it->v.z() ;
e1x = e1x + vx;
e1y = e1y + vy;
e1z = e1z + vz;
CGAL_SMS_LT_TRACE(1," vx:" << n_to_string(vx) << " vy:" << n_to_string(vy) << " vz:" << n_to_string(vz) << " e1x:"
<< n_to_string(e1x) << " e1y:" << n_to_string(e1y) << " e1z:" << n_to_string(e1z) );
}
Matrix H = LT_product(e1);
Vector c = cross_product(e1,e2);
CGAL_SMS_LT_TRACE(1
," Adding boundary preservation constraint."
<< "\n SumV:" << xyz_to_string(e1)
<< "\n SumN:" << xyz_to_string(e2)
<< "\n H:" << matrix_to_string(H)
<< "\n c:" << xyz_to_string(c)
);
Add_constraint_from_gradient(H,c);
}
}
template<class TM, class K>
void LindstromTurkCore<TM,K>::Add_volume_preservation_constraints( Triangle_data_vector const& aTriangles )
{
CGAL_SMS_LT_TRACE(1," Adding volume preservation constraint. " << aTriangles.size() << " triangles.");
Vector lSumV = NULL_VECTOR ;
FT lSumL(0) ;
for( typename Triangle_data_vector::const_iterator it = aTriangles.begin(), eit = aTriangles.end() ; it != eit ; ++it )
{
lSumV = lSumV + it->NormalV ;
lSumL = lSumL + it->NormalL ;
}
CGAL_SMS_LT_TRACE(1, " SumV:" << xyz_to_string(lSumV) << "\n SumL:" << n_to_string(lSumL) );
Add_constraint_if_alpha_compatible(lSumV,lSumL);
}
template<class TM, class K>
void LindstromTurkCore<TM,K>::Add_boundary_and_volume_optimization_constraints( Boundary_data_vector const& aBdry
, Triangle_data_vector const& aTriangles
)
{
CGAL_SMS_LT_TRACE(1," Adding boundary and volume optimization constraints. ");
Matrix H = NULL_MATRIX ;
Vector c = NULL_VECTOR ;
//
// Volume optimization
//
for( typename Triangle_data_vector::const_iterator it = aTriangles.begin(), eit = aTriangles.end() ; it != eit ; ++it )
{
Triangle_data const& lTri = *it ;
H += direct_product(lTri.NormalV,lTri.NormalV) ;
c = c - ( lTri.NormalL * lTri.NormalV ) ;
}
CGAL_SMS_LT_TRACE(2," Hv:" << matrix_to_string(H) << "\n cv:" << xyz_to_string(c) ) ;
if ( aBdry.size() > 0 )
{
//
// Boundary optimization
//
Matrix Hb = NULL_MATRIX ;
Vector cb = NULL_VECTOR ;
for ( typename Boundary_data_vector::const_iterator it = aBdry.begin() ; it != aBdry.end() ; ++ it )
{
Matrix H = LT_product(it->v);
Vector c = cross_product(it->v,it->n);
Hb += H ;
cb = cb + c ;
}
CGAL_SMS_LT_TRACE(2," Hb:" << matrix_to_string(Hb) << "\n cb:" << xyz_to_string(cb) ) ;
//
// Weighted average
//
FT lScaledBoundaryWeight = FT(9) * FT(mParams.BoundaryWeight) * squared_distance(mProfile.p0(),mProfile.p1()) ;
H *= mParams.VolumeWeight ;
c = c * mParams.VolumeWeight ;
H += lScaledBoundaryWeight * Hb ;
c = c + ( lScaledBoundaryWeight * cb ) ;
CGAL_SMS_LT_TRACE(2," H:" << matrix_to_string(H) << "\n c:" << xyz_to_string(c) ) ;
CGAL_SMS_LT_TRACE(2," VolW:" << mParams.VolumeWeight << " BdryW:" << mParams.BoundaryWeight << " ScaledBdryW:" << lScaledBoundaryWeight ) ;
}
Add_constraint_from_gradient(H,c);
}
template<class TM, class K>
void LindstromTurkCore<TM,K>::Add_shape_optimization_constraints( vertex_descriptor_vector const& aLink )
{
FT s((double)aLink.size());
Matrix H ( s,0.0,0.0
,0.0, s,0.0
,0.0,0.0, s
);
Vector c = NULL_VECTOR ;
for( typename vertex_descriptor_vector::const_iterator it = aLink.begin(), eit = aLink.end() ; it != eit ; ++it )
c = c + (ORIGIN - get_point(*it)) ;
CGAL_SMS_LT_TRACE(1," Adding shape optimization constraint. Shape vector: " << xyz_to_string(c) );
Add_constraint_from_gradient(H,c);
}
template<class TM, class K>
typename LindstromTurkCore<TM,K>::FT
LindstromTurkCore<TM,K>::Compute_boundary_cost( Vector const& v, Boundary_data_vector const& aBdry )
{
FT rCost(0);
for ( typename Boundary_data_vector::const_iterator it = aBdry.begin() ; it != aBdry.end() ; ++ it )
{
Vector u = (it->t - ORIGIN ) - v ;
Vector c = cross_product(it->v,u);
rCost += c*c;
}
return rCost / FT(4) ;
}
template<class TM, class K>
typename LindstromTurkCore<TM,K>::FT
LindstromTurkCore<TM,K>::Compute_volume_cost( Vector const& v, Triangle_data_vector const& aTriangles )
{
FT rCost(0);
for( typename Triangle_data_vector::const_iterator it = aTriangles.begin(), eit = aTriangles.end() ; it != eit ; ++it )
{
Triangle_data const& lTri = *it ;
FT lF = lTri.NormalV * v - lTri.NormalL ;
rCost += ( lF * lF ) ;
}
return rCost / FT(36) ;
}
template<class TM, class K>
typename LindstromTurkCore<TM,K>::FT
LindstromTurkCore<TM,K>::Compute_shape_cost( Point const& p, vertex_descriptor_vector const& aLink )
{
FT rCost(0);
for( typename vertex_descriptor_vector::const_iterator it = aLink.begin(), eit = aLink.end() ; it != eit ; ++it )
rCost += squared_distance(p,get_point(*it)) ;
return rCost ;
}
template<class TM, class K>
void LindstromTurkCore<TM,K>::Add_constraint_if_alpha_compatible( Vector const& Ai, FT const& bi )
{
CGAL_SMS_LT_TRACE(3," Adding new constraints if alpha-compatible.\n Ai: " << xyz_to_string(Ai) << "\n bi:" << n_to_string(bi) << ")" );
FT slai = Ai*Ai ;
CGAL_SMS_LT_TRACE(3,"\n slai: " << n_to_string(slai) << ")" );
if ( is_finite(slai) )
{
FT l = CGAL_NTS sqrt( slai ) ;
CGAL_SMS_LT_TRACE(3," l: " << n_to_string(l) );
if ( !CGAL_NTS is_zero(l) )
{
Vector Ain = Ai / l ;
FT bin = bi / l ;
CGAL_SMS_LT_TRACE(3," Ain: " << xyz_to_string(Ain) << " bin:" << n_to_string(bin) );
bool lAddIt = true ;
if ( mConstraints_n == 1 )
{
FT d01 = mConstraints_A.r0() * Ai ;
FT sla0 = mConstraints_A.r0() * mConstraints_A.r0() ;
FT sd01 = d01 * d01 ;
FT max = sla0 * slai * mSquared_cos_alpha ;
CGAL_SMS_LT_TRACE(3," Second constraint. d01: " << n_to_string(d01) << " sla0:" << n_to_string(sla0) << " sd01:" << n_to_string(sd01) << " max:" << n_to_string(max) );
if ( sd01 > max )
lAddIt = false ;
}
else if ( mConstraints_n == 2 )
{
Vector N = cross_product(mConstraints_A.r0(),mConstraints_A.r1());
FT dc012 = N * Ai ;
FT slc01 = N * N ;
FT sdc012 = dc012 * dc012;
FT min = slc01 * slai * mSquared_sin_alpha ;
CGAL_SMS_LT_TRACE(3," Third constraint. N: " << xyz_to_string(N) << " dc012:" << n_to_string(dc012) << " slc01:" << n_to_string(slc01)
<< " sdc012:" << n_to_string(sdc012) << " min:" << n_to_string(min) );
if ( sdc012 <= min )
lAddIt = false ;
}
if ( lAddIt )
{
switch ( mConstraints_n )
{
case 0 :
mConstraints_A.r0() = Ain ;
mConstraints_b = Vector(bin,mConstraints_b.y(),mConstraints_b.z());
break ;
case 1 :
mConstraints_A.r1() = Ain ;
mConstraints_b = Vector(mConstraints_b.x(),bin,mConstraints_b.z());
break ;
case 2 :
mConstraints_A.r2() = Ain ;
mConstraints_b = Vector(mConstraints_b.x(),mConstraints_b.y(),bin);
break ;
}
CGAL_SMS_LT_TRACE(3," Accepting # " << mConstraints_n << " A:" << matrix_to_string(mConstraints_A) << " b:" << xyz_to_string(mConstraints_b) ) ;
++ mConstraints_n ;
}
else
{
CGAL_SMS_LT_TRACE(3," INCOMPATIBLE. Discarded" ) ;
}
}
else
{
CGAL_SMS_LT_TRACE(3," l is ZERO. Discarded" );
}
}
else
{
CGAL_SMS_LT_TRACE(3," OVERFLOW. Discarded" ) ;
}
}
template<class V>
int index_of_max_component ( V const& v )
{
typedef typename Kernel_traits<V>::Kernel::FT FT ;
int i = 0 ;
FT max = v.x();
if ( max < v.y() )
{
max = v.y();
i = 1 ;
}
if ( max < v.z() )
{
max = v.z();
i = 2 ;
}
return i ;
}
template<class TM, class K>
void LindstromTurkCore<TM,K>::Add_constraint_from_gradient ( Matrix const& H, Vector const& c )
{
CGAL_SMS_LT_TRACE(3," Adding constraint from gradient. Current n=" << mConstraints_n ) ;
CGAL_precondition(mConstraints_n >= 0 && mConstraints_n<=2 );
switch(mConstraints_n)
{
case 0 :
Add_constraint_if_alpha_compatible(H.r0(),-c.x());
Add_constraint_if_alpha_compatible(H.r1(),-c.y());
Add_constraint_if_alpha_compatible(H.r2(),-c.z());
break;
case 1 :
{
Vector const& A0 = mConstraints_A.r0();
CGAL_assertion( A0 != NULL_VECTOR ) ;
Vector AbsA0( CGAL::abs(A0.x())
, CGAL::abs(A0.y())
, CGAL::abs(A0.z())
);
Vector Q0;
switch ( index_of_max_component(AbsA0) )
{
// Since A0 is guaranteed to be non-zero, the denominators here are known to be non-zero too.
case 0: Q0 = Vector(- A0.z()/A0.x(),0 ,1 ); break;
case 1: Q0 = Vector(0 ,- A0.z()/A0.y(),1 ); break;
case 2: Q0 = Vector(1 ,0 ,- A0.x()/A0.z()); break;
default : Q0 = NULL_VECTOR ; // This should never happen!
}
CGAL_SMS_LT_TRACE(3," Q0:" << xyz_to_string(Q0) ) ;
CGAL_assertion( Q0 != NULL_VECTOR ) ;
Vector Q1 = cross_product(A0,Q0);
Vector A1 = H * Q0 ;
Vector A2 = H * Q1 ;
FT b1 = - ( Q0 * c ) ;
FT b2 = - ( Q1 * c ) ;
CGAL_SMS_LT_TRACE(3," Q1:" << xyz_to_string(Q1) << "\n A1: " << xyz_to_string(A1) << "\n A2:" << xyz_to_string(A2)
<< "\n b1:" << n_to_string(b1) << "\n b2:" << n_to_string(b2) ) ;
Add_constraint_if_alpha_compatible(A1,b1);
Add_constraint_if_alpha_compatible(A2,b2);
}
break ;
case 2:
{
Vector Q = cross_product(mConstraints_A.r0(),mConstraints_A.r1());
Vector A2 = H * Q ;
FT b2 = - ( Q * c ) ;
CGAL_SMS_LT_TRACE(3," Q:" << xyz_to_string(Q) << "\n A2: " << xyz_to_string(A2) << "\n b2:" << n_to_string(b2) ) ;
Add_constraint_if_alpha_compatible(A2,b2);
}
break ;
}
}
} // namespace Surface_mesh_simplification
} //namespace CGAL
#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_EDGE_COLLAPSE_DETAIL_LINDSTROMTURK_CORE_IMPL_H //
// EOF //

33
Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_length_cost.h
View File

@ -22,51 +22,28 @@
#include <CGAL/license/Surface_mesh_simplification.h>
#include <CGAL/Surface_mesh_simplification/Detail/Common.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_profile.h>
namespace CGAL {
namespace Surface_mesh_simplification {
namespace Surface_mesh_simplification
{
//
// Edge-length cost: the squared length of the collapsing edge
//
template<class TM>
template<class TM>
class Edge_length_cost
{
public:
/*
typedef TM_ TM ;
typedef Edge_profile<TM> Profile ;
typedef typename Profile::Point Point;
typedef typename Kernel_traits<Point>::Kernel Kernel ;
typedef typename Kernel::FT FT ;
typedef optional<FT> result_type ;
*/
public:
Edge_length_cost()
{}
Edge_length_cost() {}
template <typename Profile, typename T>
optional<typename Profile::FT> operator()( Profile const& aProfile, T const& /*aPlacement*/ ) const
optional<typename Profile::FT> operator()(const Profile& aProfile, const T& /*aPlacement*/) const
{
typedef optional<typename Profile::FT> result_type;
return result_type(squared_distance(aProfile.p0(),aProfile.p1()));
return result_type(squared_distance(aProfile.p0(), aProfile.p1()));
}
};
} // namespace Surface_mesh_simplification
} //namespace CGAL
#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_EDGE_COLLAPSE_EDGE_LENGHT_COST_H
// EOF //

25
Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_length_stop_predicate.h
View File

@ -24,34 +24,29 @@
#include <CGAL/squared_distance_3.h>
namespace CGAL {
namespace Surface_mesh_simplification
{
namespace Surface_mesh_simplification {
template <class FT>
class Edge_length_stop_predicate
{
FT m_edge_sq_length_threshold;
public:
Edge_length_stop_predicate( double edge_length_threshold )
: m_edge_sq_length_threshold(edge_length_threshold*edge_length_threshold)
Edge_length_stop_predicate(double edge_length_threshold)
: m_edge_sq_length_threshold(edge_length_threshold * edge_length_threshold)
{}
template <typename F, typename Profile>
bool operator()( F const&
, Profile const& profile
, std::size_t
, std::size_t
) const
bool operator()(const F&,
const Profile& profile,
std::size_t,
std::size_t) const
{
return CGAL::squared_distance(profile.p0(), profile.p1()) >
m_edge_sq_length_threshold;
return CGAL::squared_distance(profile.p0(), profile.p1()) > m_edge_sq_length_threshold;
}
};
} // namespace Surface_mesh_simplification
} //namespace CGAL
#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_EDGE_COLLAPSE_EDGE_LENGTH_STOP_PREDICATE_H //
#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_EDGE_COLLAPSE_EDGE_LENGTH_STOP_PREDICATE_H

353
Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_profile.h
View File

@ -18,160 +18,140 @@
// Author(s) : Fernando Cacciola <fernando.cacciola@geometryfactory.com>
//
#ifndef CGAL_SURFACE_MESH_SIMPLIFICATION_DETAIL_EDGE_PROFILE_H
#define CGAL_SURFACE_MESH_SIMPLIFICATION_DETAIL_EDGE_PROFILE_H 1
#define CGAL_SURFACE_MESH_SIMPLIFICATION_DETAIL_EDGE_PROFILE_H
#include <CGAL/license/Surface_mesh_simplification.h>
#include <CGAL/Surface_mesh_simplification/Detail/Common.h>
#include <vector>
#include <set>
#include <CGAL/Surface_mesh_simplification/Detail/Common.h>
namespace CGAL {
namespace Surface_mesh_simplification {
namespace Surface_mesh_simplification
{
template<class TM_, class VertexPointMap_ = typename boost::property_map<TM_, CGAL::vertex_point_t>::type>
template<class TM_,
class VertexPointMap_ = typename boost::property_map<TM_, CGAL::vertex_point_t>::type>
class Edge_profile
{
public:
typedef TM_ TM ;
typedef TM_ TM;
typedef VertexPointMap_ VertexPointMap;
typedef boost::graph_traits<TM> GraphTraits ;
typedef boost::graph_traits<TM> GraphTraits;
typedef typename GraphTraits::vertex_descriptor vertex_descriptor ;
typedef typename GraphTraits::face_descriptor face_descriptor ;
typedef typename GraphTraits::halfedge_descriptor halfedge_descriptor ;
typedef typename GraphTraits::vertex_descriptor vertex_descriptor;
typedef typename GraphTraits::halfedge_descriptor halfedge_descriptor;
typedef typename GraphTraits::face_descriptor face_descriptor;
typedef std::vector<vertex_descriptor> vertex_descriptor_vector;
typedef std::vector<halfedge_descriptor> halfedge_descriptor_vector;
//typedef typename boost::property_map<TM, CGAL::vertex_point_t>::type Vertex_point_pmap;
typedef typename boost::property_traits<VertexPointMap>::value_type Point;
typedef typename Kernel_traits<Point>::Kernel Kernel;
typedef typename Kernel::FT FT;
public:
struct Triangle
{
Triangle() {}
Triangle( vertex_descriptor const& a_v0
, vertex_descriptor const& a_v1
, vertex_descriptor const& a_v2
)
Triangle(const vertex_descriptor a_v0,
const vertex_descriptor a_v1,
const vertex_descriptor a_v2)
: v0(a_v0), v1(a_v1), v2(a_v2)
{
CGAL_SURF_SIMPL_TEST_assertion( handle_assigned(v0) && handle_assigned(v1) && handle_assigned(v2) ) ;
CGAL_SURF_SIMPL_TEST_assertion( v0 != v1 && v1 != v2 && v2 != v0 ) ;
CGAL_SURF_SIMPL_TEST_assertion(handle_assigned(v0) && handle_assigned(v1) && handle_assigned(v2));
CGAL_SURF_SIMPL_TEST_assertion(v0 != v1 && v1 != v2 && v2 != v0);
}
vertex_descriptor v0 ;
vertex_descriptor v1 ;
vertex_descriptor v2 ;
} ;
vertex_descriptor v0;
vertex_descriptor v1;
vertex_descriptor v2;
};
typedef std::vector<vertex_descriptor> vertex_descriptor_vector ;
typedef std::vector<halfedge_descriptor> halfedge_descriptor_vector ;
typedef std::vector<Triangle> Triangle_vector ;
typedef std::vector<Triangle> Triangle_vector;
public :
template<class VertexIdxMap
,class EdgeIdxMap
>
Edge_profile ( halfedge_descriptor const& aV0V1
, TM& aSurface
, VertexIdxMap const& aVertex_index_map
, VertexPointMap const& aVertex_point_map
, EdgeIdxMap const& aEdge_index_map
, bool has_border
) ;
template<class VertexIdxMap,
class EdgeIdxMap>
Edge_profile(const halfedge_descriptor& aV0V1,
TM& aSurface,
const VertexIdxMap& aVertex_index_map,
const VertexPointMap& aVertex_point_map,
const EdgeIdxMap& aEdge_index_map,
bool has_border);
public :
const halfedge_descriptor v0_v1() const { return mV0V1; }
const halfedge_descriptor v1_v0() const { return mV1V0; }
halfedge_descriptor const& v0_v1() const { return mV0V1; }
halfedge_descriptor const& v1_v0() const { return mV1V0; }
vertex_descriptor const& v0() const { return mV0; }
vertex_descriptor const& v1() const { return mV1; }
const vertex_descriptor v0() const { return mV0; }
const vertex_descriptor v1() const { return mV1; }
// These are null if v0v1 is a border (thius there is no face to its left)
vertex_descriptor const& vL() const { return mVL; }
halfedge_descriptor const& v1_vL() const { return mV1VL; }
halfedge_descriptor const& vL_v0() const { return mVLV0; }
const vertex_descriptor vL() const { return mVL; }
const halfedge_descriptor v1_vL() const { return mV1VL; }
const halfedge_descriptor vL_v0() const { return mVLV0; }
// These are null if v1v0 is a border (thius there is no face to its left)
vertex_descriptor const& vR() const { return mVR; }
halfedge_descriptor const& v0_vR() const { return mV0VR; }
halfedge_descriptor const& vR_v1() const { return mVRV1; }
const vertex_descriptor vR() const { return mVR; }
const halfedge_descriptor v0_vR() const { return mV0VR; }
const halfedge_descriptor vR_v1() const { return mVRV1; }
Triangle_vector const& triangles() const {
const Triangle_vector& triangles() const
{
if(mTriangles.empty())
const_cast<Edge_profile*>(this)->extract_triangles_and_link();
if(mTriangles.empty()){
const_cast<Edge_profile*>(this)->Extract_triangles_and_link();
}
CGAL_HISTOGRAM_PROFILER("triangles.size()", mTriangles.size());
return mTriangles ; }
return mTriangles;
}
// The cycle of vertices around the edge
vertex_descriptor_vector const& link() const {
const vertex_descriptor_vector& link() const
{
CGAL_PROFILER("link calls");
if(mLink.empty()){
if(mLink.empty())
const_cast<Edge_profile*>(this)->extract_triangles_and_link();
const_cast<Edge_profile*>(this)->Extract_triangles_and_link();
return mLink;
}
return mLink ; }
halfedge_descriptor_vector const& border_edges() const {
return mBorderEdges ;
}
TM& surface() const { return *mSurface ; }
TM& surface_mesh() const { return *mSurface ; }
VertexPointMap vertex_point_map() const { return mvpm ; }
const halfedge_descriptor_vector& border_edges() const { return mBorderEdges; }
TM& surface() const { return *mSurface; }
TM& surface_mesh() const { return *mSurface; }
VertexPointMap vertex_point_map() const { return mvpm; }
public :
const Point& p0() const { return mP0; }
const Point& p1() const { return mP1; }
Point const& p0() const { return mP0; }
Point const& p1() const { return mP1; }
bool is_v0_v1_a_border() const { return mIsBorderV0V1; }
bool is_v1_v0_a_border() const { return mIsBorderV1V0; }
bool is_v0_v1_a_border() const { return mIsBorderV0V1 ; }
bool is_v1_v0_a_border() const { return mIsBorderV1V0 ; }
bool left_face_exists () const { return !mIsBorderV0V1 ; }
bool right_face_exists() const { return !mIsBorderV1V0 ; }
bool left_face_exists () const { return !mIsBorderV0V1; }
bool right_face_exists() const { return !mIsBorderV1V0; }
private:
// typedef typename GraphTraits::in_edge_iterator in_edge_iterator ;
bool is_border(halfedge_descriptor e) const
{
bool is_border(halfedge_descriptor e) const {
return face(e,*mSurface) == boost::graph_traits<TM>::null_face();
}
void Extract_borders() ;
void Extract_triangles_and_link() ;
void extract_borders();
void extract_triangles_and_link();
private:
halfedge_descriptor mV0V1;
halfedge_descriptor mV1V0;
bool mIsBorderV0V1 ;
bool mIsBorderV1V0 ;
bool mIsBorderV0V1;
bool mIsBorderV1V0;
vertex_descriptor mV0;
vertex_descriptor mV1;
Point mP0 ;
Point mP1 ;
Point mP0;
Point mP1;
vertex_descriptor mVL;
vertex_descriptor mVR;
@ -181,20 +161,197 @@ private:
halfedge_descriptor mV0VR;
halfedge_descriptor mVRV1;
vertex_descriptor_vector mLink ;
halfedge_descriptor_vector mBorderEdges ;
Triangle_vector mTriangles ;
vertex_descriptor_vector mLink;
halfedge_descriptor_vector mBorderEdges;
Triangle_vector mTriangles;
TM* mSurface ;
TM* mSurface;
VertexPointMap mvpm;
} ;
};
template<class TM, class VertexPointMap>
template<class VertexIdxMap, class EdgeIdxMap>
Edge_profile<TM, VertexPointMap>::
Edge_profile(const halfedge_descriptor& aV0V1,
TM& aSurface,
const VertexIdxMap&,
const VertexPointMap& aVertex_point_map,
const EdgeIdxMap&,
bool has_border)
:
mV0V1(aV0V1),
mSurface(boost::addressof(aSurface)),
mvpm(aVertex_point_map)
{
CGAL_PROFILER("Edge_profile constructor calls");
mLink.reserve(12);
mTriangles.reserve(16);
mV1V0 = opposite(v0_v1(), surface_mesh());
mV0 = source(v0_v1(), surface_mesh());
mV1 = target(v0_v1(), surface_mesh());
CGAL_assertion(mV0 != mV1);
mP0 = get(vertex_point_map(),mV0);
mP1 = get(vertex_point_map(),mV1);
mIsBorderV0V1 = is_border(v0_v1());
mIsBorderV1V0 = is_border(v1_v0());
if(left_face_exists())
{
CGAL_SURF_SIMPL_TEST_assertion(! is_border(mV0V1));
mVLV0 = prev(v0_v1(), surface_mesh());
mV1VL = next(v0_v1(), surface_mesh());
mVL = target(v1_vL(), surface_mesh());
CGAL_SURF_SIMPL_TEST_assertion(mV0 != mVL);
CGAL_SURF_SIMPL_TEST_assertion(mVL == source(vL_v0(), surface_mesh()));
}
else
{
CGAL_SURF_SIMPL_TEST_assertion(is_border(mV0V1));
}
if(right_face_exists())
{
CGAL_SURF_SIMPL_TEST_assertion(! is_border(mV1V0));
mV0VR = next(v1_v0(), surface_mesh());
mVRV1 = prev(v1_v0(), surface_mesh());
mVR = target(v0_vR(), surface_mesh());
CGAL_SURF_SIMPL_TEST_assertion(mV0 != mVR);
CGAL_SURF_SIMPL_TEST_assertion(mVR == source(vR_v1(), surface_mesh()));
}
else
{
CGAL_SURF_SIMPL_TEST_assertion(is_border(mV1V0));
}
if(has_border)
extract_borders();
}
template<class TM, class VertexPointMap>
void
Edge_profile<TM,VertexPointMap>::
extract_borders()
{
halfedge_descriptor e = mV0V1;
halfedge_descriptor oe = opposite(e, surface_mesh());
bool b;
if((b = is_border(e)) || is_border(oe))
mBorderEdges.push_back(b?e:oe);
e = next(oe, surface_mesh());
oe = opposite(e, surface_mesh());
while(e != mV0V1)
{
if((b = is_border(e)) || is_border(oe))
mBorderEdges.push_back(b?e:oe);
e = next(oe, surface_mesh());
oe = opposite(e, surface_mesh());
}
e = opposite(next(e, surface_mesh()), surface_mesh());
oe = opposite(e, surface_mesh());
while(e != mV0V1)
{
if((b = is_border(e)) || is_border(oe))
mBorderEdges.push_back(b?e:oe);
e = opposite(next(e, surface_mesh()), surface_mesh());
oe = opposite(e, surface_mesh());
}
}
// Extract all triangles (its normals) and vertices (the link) around the collapsing edge p_q
template<class TM, class VertexPointMap>
void
Edge_profile<TM,VertexPointMap>::
extract_triangles_and_link()
{
#ifdef CGAL_SMS_EDGE_PROFILE_ALWAYS_NEED_UNIQUE_VERTEX_IN_LINK
std::set<vertex_descriptor> vertex_already_inserted;
#endif
// look at the two faces or holes adjacent to edge (v0,v1)
// and at the opposite vertex if it exists
halfedge_descriptor endleft = next(v1_v0(), surface_mesh());
halfedge_descriptor endright = next(v0_v1(), surface_mesh());
if(left_face_exists())
mTriangles.push_back(Triangle(v0(), v1(), vL()));
if(right_face_exists())
mTriangles.push_back(Triangle(v1(), v0(), vR()));
// counterclockwise around v0
halfedge_descriptor e02 = opposite(prev(v0_v1(), surface_mesh()), surface_mesh());
vertex_descriptor v = target(e02, surface_mesh()), v2 = v;
while(e02 != endleft)
{
#ifdef CGAL_SMS_EDGE_PROFILE_ALWAYS_NEED_UNIQUE_VERTEX_IN_LINK
if(vertex_already_inserted.insert(v2).second)
#endif
mLink.push_back(v2);
bool is_b = is_border(e02);
e02 = opposite(prev(e02, surface_mesh()), surface_mesh());
v = target(e02, surface_mesh());
if(!is_b)
mTriangles.push_back(Triangle(v, v0(), v2));
v2 = v;
}
e02 = opposite(prev(v1_v0(), surface_mesh()), surface_mesh());
if(target(e02, surface_mesh()) != v) // add the vertex if it is not added in the following loop
{
#ifdef CGAL_SMS_EDGE_PROFILE_ALWAYS_NEED_UNIQUE_VERTEX_IN_LINK
if(vertex_already_inserted.insert(v).second)
#endif
mLink.push_back(v);
}
// counterclockwise around v1
v2 = target(e02, surface_mesh());
v = v2;
while(e02 != endright)
{
#ifdef CGAL_SMS_EDGE_PROFILE_ALWAYS_NEED_UNIQUE_VERTEX_IN_LINK
if(vertex_already_inserted.insert(v2).second)
#endif
mLink.push_back(v2);
bool is_b = is_border(e02);
e02 = opposite(prev(e02, surface_mesh()), surface_mesh());
v = target(e02, surface_mesh());
if(!is_b)
mTriangles.push_back(Triangle(v,v1(),v2));
v2 = v;
}
if(mLink.empty() || //handle link of an isolated triangle
target(opposite(prev(v0_v1(), surface_mesh()), surface_mesh()), surface_mesh())!=v)
{
#ifdef CGAL_SMS_EDGE_PROFILE_ALWAYS_NEED_UNIQUE_VERTEX_IN_LINK
if(vertex_already_inserted.insert(v).second)
#endif
mLink.push_back(v);
}
CGAL_assertion(!mLink.empty());
}
} // namespace Surface_mesh_simplification
} //namespace CGAL
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_profile_impl.h>
#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_DETAIL_EDGE_PROFILE_H
// EOF //

214
Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_profile_impl.h
View File

@ -1,214 +0,0 @@
// Copyright (c) 2006 GeometryFactory (France). All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
// You can redistribute it and/or modify it under the terms of the GNU
// General Public License as published by the Free Software Foundation,
// either version 3 of the License, or (at your option) any later version.
//
// Licensees holding a valid Surface_mesh_simplification license may use this file in
// accordance with the Surface_mesh_simplification license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0+
//
// Author(s) : Fernando Cacciola <fernando.cacciola@geometryfactory.com>
//
#ifndef CGAL_SURFACE_MESH_SIMPLIFICATION_DETAIL_EDGE_PROFILE_IMPL_H
#define CGAL_SURFACE_MESH_SIMPLIFICATION_DETAIL_EDGE_PROFILE_IMPL_H 1
#include <CGAL/license/Surface_mesh_simplification.h>
namespace CGAL {
namespace Surface_mesh_simplification
{
template<class TM, class VertexPointMap>
template<class VertexIdxMap
,class EdgeIdxMap
>
Edge_profile<TM,VertexPointMap>::Edge_profile ( halfedge_descriptor const& aV0V1
, TM& aSurface
, VertexIdxMap const&
, VertexPointMap const& aVertex_point_map
, EdgeIdxMap const&
, bool has_border
)
:
mV0V1(aV0V1)
,mSurface(boost::addressof(aSurface))
, mvpm(aVertex_point_map)
{
CGAL_PROFILER("Edge_profile constructor calls");
mLink.reserve(12);
mTriangles.reserve(16);
mV1V0 = opposite(v0_v1(),surface_mesh());
mV0 = source(v0_v1(),surface_mesh());
mV1 = target(v0_v1(),surface_mesh());
CGAL_assertion( mV0 != mV1 );
mP0 = get(vertex_point_map(),mV0);
mP1 = get(vertex_point_map(),mV1);
mIsBorderV0V1 = is_border(v0_v1());
mIsBorderV1V0 = is_border(v1_v0());
if ( left_face_exists() )
{
CGAL_SURF_SIMPL_TEST_assertion( ! is_border(mV0V1) ) ;
mVLV0 = prev(v0_v1(),surface_mesh());
mV1VL = next(v0_v1(),surface_mesh());
mVL = target(v1_vL(),surface_mesh());
CGAL_SURF_SIMPL_TEST_assertion( mV0 != mVL );
CGAL_SURF_SIMPL_TEST_assertion( mVL == source(vL_v0(),surface_mesh()) );
}
else
{
CGAL_SURF_SIMPL_TEST_assertion( is_border(mV0V1) ) ;
}
if ( right_face_exists() )
{
CGAL_SURF_SIMPL_TEST_assertion( ! is_border(mV1V0) ) ;
mV0VR = next(v1_v0(),surface_mesh());
mVRV1 = prev(v1_v0(),surface_mesh());
mVR = target(v0_vR(),surface_mesh());
CGAL_SURF_SIMPL_TEST_assertion( mV0 != mVR );
CGAL_SURF_SIMPL_TEST_assertion( mVR == source(vR_v1(),surface_mesh()) );
}
else
{
CGAL_SURF_SIMPL_TEST_assertion( is_border(mV1V0) ) ;
}
if(has_border){
Extract_borders();
}
}
template<class TM, class VertexPointMap>
void Edge_profile<TM,VertexPointMap>::Extract_borders()
{
halfedge_descriptor e = mV0V1;
halfedge_descriptor oe = opposite(e, surface_mesh());
bool b;
if((b = is_border(e)) || is_border(oe)){
mBorderEdges.push_back(b?e:oe);
}
e = next(oe,surface_mesh());
oe = opposite(e,surface_mesh());
while(e != mV0V1){
if((b = is_border(e)) || is_border(oe)){
mBorderEdges.push_back(b?e:oe);
}
e = next(oe,surface_mesh());
oe = opposite(e,surface_mesh());
}
e = opposite(next(e,surface_mesh()),surface_mesh());
oe = opposite(e,surface_mesh());
while(e != mV0V1){
if((b = is_border(e)) || is_border(oe)){
mBorderEdges.push_back(b?e:oe);
}
e = opposite(next(e,surface_mesh()),surface_mesh());
oe = opposite(e,surface_mesh());
}
}
// Extract all triangles (its normals) and vertices (the link) around the collapsing edge p_q
//
template<class TM, class VertexPointMap>
void Edge_profile<TM,VertexPointMap>::Extract_triangles_and_link()
{
#ifdef CGAL_SMS_EDGE_PROFILE_ALWAYS_NEED_UNIQUE_VERTEX_IN_LINK
std::set<vertex_descriptor> vertex_already_inserted;
#endif
// look at the two faces or holes adjacent to edge (v0,v1)
// and at the opposite vertex if it exists
halfedge_descriptor endleft = next(v1_v0(), surface_mesh());
halfedge_descriptor endright = next(v0_v1(), surface_mesh());
if( left_face_exists() )
mTriangles.push_back(Triangle(v0(),v1(),vL()) ) ;
if( right_face_exists() )
mTriangles.push_back(Triangle(v1(),v0(),vR()) ) ;
// counterclockwise around v0
halfedge_descriptor e02 = opposite(prev(v0_v1(),surface_mesh()), surface_mesh());
vertex_descriptor v=target(e02,surface_mesh()), v2=v;
while(e02 != endleft) {
#ifdef CGAL_SMS_EDGE_PROFILE_ALWAYS_NEED_UNIQUE_VERTEX_IN_LINK
if (vertex_already_inserted.insert(v2).second)
#endif
mLink.push_back(v2);
bool is_b = is_border(e02);
e02 = opposite(prev(e02,surface_mesh()), surface_mesh());
v = target(e02,surface_mesh());
if( !is_b )
mTriangles.push_back(Triangle(v,v0(),v2) ) ;
v2 = v;
}
e02 = opposite(prev(v1_v0(),surface_mesh()), surface_mesh());
if(target(e02, surface_mesh())!=v ) // add the vertex if it is not added in the following loop
{
#ifdef CGAL_SMS_EDGE_PROFILE_ALWAYS_NEED_UNIQUE_VERTEX_IN_LINK
if (vertex_already_inserted.insert(v).second)
#endif
mLink.push_back(v);
}
// counterclockwise around v1
v2 = target(e02,surface_mesh());
v = v2;
while(e02 != endright) {
#ifdef CGAL_SMS_EDGE_PROFILE_ALWAYS_NEED_UNIQUE_VERTEX_IN_LINK
if (vertex_already_inserted.insert(v2).second)
#endif
mLink.push_back(v2);
bool is_b = is_border(e02);
e02 = opposite(prev(e02,surface_mesh()), surface_mesh());
v = target(e02,surface_mesh());
if( !is_b ){
mTriangles.push_back(Triangle(v,v1(),v2) ) ;
}
v2 = v;
}
if(mLink.empty() || //handle link of an isolated triangle
target(opposite(prev(v0_v1(),surface_mesh()), surface_mesh()), surface_mesh())!=v)
{
#ifdef CGAL_SMS_EDGE_PROFILE_ALWAYS_NEED_UNIQUE_VERTEX_IN_LINK
if (vertex_already_inserted.insert(v).second)
#endif
mLink.push_back(v);
}
CGAL_assertion(!mLink.empty());
}
} // namespace Surface_mesh_simplification
} //namespace CGAL
#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_DETAIL_EDGE_PROFILE_IMPL_H
// EOF //

7
Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/LindstromTurk.h
View File

@ -18,15 +18,12 @@
// Author(s) : Fernando Cacciola <fernando.cacciola@geometryfactory.com>
//
#ifndef CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_EDGE_COLLAPSE_LINDSTROMTURK_H
#define CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_EDGE_COLLAPSE_LINDSTROMTURK_H 1
#define CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_EDGE_COLLAPSE_LINDSTROMTURK_H
#include <CGAL/license/Surface_mesh_simplification.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/LindstromTurk_params.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/LindstromTurk_cost.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/LindstromTurk_placement.h>
#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_EDGE_COLLAPSE_LINDSTROMTURK_H //
// EOF //
#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_EDGE_COLLAPSE_LINDSTROMTURK_H

38
Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/LindstromTurk_cost.h
View File

@ -18,53 +18,39 @@
// Author(s) : Fernando Cacciola <fernando.cacciola@geometryfactory.com>
//
#ifndef CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_EDGE_COLLAPSE_LINDSTROMTURK_COST_H
#define CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_EDGE_COLLAPSE_LINDSTROMTURK_COST_H 1
#define CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_EDGE_COLLAPSE_LINDSTROMTURK_COST_H
#include <CGAL/license/Surface_mesh_simplification.h>
#include <CGAL/Surface_mesh_simplification/Detail/Common.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Detail/Lindstrom_Turk_core.h>
namespace CGAL {
namespace Surface_mesh_simplification {
namespace Surface_mesh_simplification
{
template<class TM_>
template<class TM_>
class LindstromTurk_cost
{
public:
typedef TM_ TM;
typedef TM_ TM ;
/*
typedef Edge_profile<TM> Profile ;
typedef typename Traits::Point_3 Point;
typedef typename Traits::FT FT ;
typedef optional<FT> result_type ;
*/
public:
LindstromTurk_cost( LindstromTurk_params const& aParams = LindstromTurk_params() ) : mParams(aParams) {}
LindstromTurk_cost(LindstromTurk_params const& aParams = LindstromTurk_params())
: mParams(aParams)
{}
template <typename Profile>
optional<typename Profile::FT>
operator()( Profile const& aProfile, optional<typename Profile::Point> const& aPlacement ) const
operator()(const Profile& aProfile,
const optional<typename Profile::Point>& aPlacement) const
{
return LindstromTurkCore<TM,Profile>(mParams,aProfile).compute_cost(aPlacement) ;
return LindstromTurkCore<TM,Profile>(mParams,aProfile).compute_cost(aPlacement);
}
private:
LindstromTurk_params mParams ;
LindstromTurk_params mParams;
};
} // namespace Surface_mesh_simplification
} //namespace CGAL
#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_EDGE_COLLAPSE_LINDSTROMTURK_COST_H //
// EOF //
#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_EDGE_COLLAPSE_LINDSTROMTURK_COST_H

30
Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/LindstromTurk_params.h
View File

@ -22,37 +22,31 @@
#include <CGAL/license/Surface_mesh_simplification.h>
namespace CGAL {
namespace Surface_mesh_simplification
{
namespace Surface_mesh_simplification {
struct LindstromTurk_params
{
LindstromTurk_params()
:
VolumeWeight (0.5)
,BoundaryWeight(0.5)
,ShapeWeight (0)
VolumeWeight(0.5),
BoundaryWeight(0.5),
ShapeWeight(0)
{}
LindstromTurk_params( double aVolumeWeight, double aBoundaryWeight, double aShapeWeight )
LindstromTurk_params(double aVolumeWeight, double aBoundaryWeight, double aShapeWeight)
:
VolumeWeight (aVolumeWeight)
,BoundaryWeight(aBoundaryWeight)
,ShapeWeight (aShapeWeight)
VolumeWeight(aVolumeWeight),
BoundaryWeight(aBoundaryWeight),
ShapeWeight(aShapeWeight)
{}
double VolumeWeight ;
double BoundaryWeight ;
double ShapeWeight ;
double VolumeWeight;
double BoundaryWeight;
double ShapeWeight;
};
} // namespace Surface_mesh_simplification
} //namespace CGAL
} // namespace CGAL
#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_EDGE_COLLAPSE_LINDSTROMTURK_PARAMS_H
// EOF //

31
Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/LindstromTurk_placement.h
View File

@ -22,44 +22,33 @@
#include <CGAL/license/Surface_mesh_simplification.h>
#include <CGAL/Surface_mesh_simplification/Detail/Common.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Detail/Lindstrom_Turk_core.h>
namespace CGAL {
namespace Surface_mesh_simplification {
namespace Surface_mesh_simplification
{
template<class TM_>
template<class TM_>
class LindstromTurk_placement
{
public:
typedef TM_ TM;
typedef TM_ TM ;
public:
LindstromTurk_placement( LindstromTurk_params const& aParams = LindstromTurk_params() ) : mParams(aParams) {}
LindstromTurk_placement(LindstromTurk_params const& aParams = LindstromTurk_params())
: mParams(aParams)
{}
template <typename Profile>
optional<typename Profile::Point>
operator()( Profile const& aProfile) const
optional<typename Profile::Point> operator()(const Profile& aProfile) const
{
return LindstromTurkCore<TM,Profile>(mParams,aProfile).compute_placement() ;
return LindstromTurkCore<TM,Profile>(mParams,aProfile).compute_placement();
}
private:
LindstromTurk_params mParams ;
LindstromTurk_params mParams;
};
} // namespace Surface_mesh_simplification
} //namespace CGAL
#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_EDGE_COLLAPSE_LINDSTROMTURK_PLACEMENT_H //
// EOF //
#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_EDGE_COLLAPSE_LINDSTROMTURK_PLACEMENT_H

3
Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Midpoint_and_length.h
View File

@ -22,10 +22,7 @@
#include <CGAL/license/Surface_mesh_simplification.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_length_cost.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Midpoint_placement.h>
#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_EDGE_COLLAPSE_MIDPOINT_AND_LENGTH_H
// EOF //

2
Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Midpoint_placement.h
View File

@ -42,7 +42,7 @@ public:
{}
template <typename Profile>
optional<typename Profile::Point> operator()( Profile const& aProfile ) const
optional<typename Profile::Point> operator()(const Profile& aProfile) const
{
return optional<typename Profile::Point> (midpoint(aProfile.p0(),aProfile.p1()));
}

188
Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/edge_collapse.h
View File

@ -22,7 +22,6 @@
#include <CGAL/license/Surface_mesh_simplification.h>
#include <CGAL/boost/graph/properties.h>
#include <CGAL/boost/graph/named_function_params.h>
@ -31,143 +30,118 @@
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/LindstromTurk.h>
namespace CGAL {
namespace Surface_mesh_simplification {
namespace Surface_mesh_simplification
{
template<class TM
,class ShouldStop
,class VertexIndexMap
,class VertexPointMap
,class EdgeIndexMap
,class EdgeIsConstrainedMap
,class GetCost
,class GetPlacement
,class Visitor
>
int edge_collapse ( TM& aSurface
, ShouldStop const& aShould_stop
template<class TM,
class ShouldStop,
class VertexIndexMap,
class VertexPointMap,
class EdgeIndexMap,
class EdgeIsConstrainedMap,
class GetCost,
class GetPlacement,
class Visitor>
int edge_collapse(TM& aSurface,
const ShouldStop& aShould_stop,
// optional mesh information policies
, VertexIndexMap const& aVertex_index_map // defaults to get(vertex_index,aSurface)
, VertexPointMap const& aVertex_point_map // defaults to get(vertex_point,aSurface)
, EdgeIndexMap const& aEdge_index_map // defaults to get(edge_index,aSurface)
, EdgeIsConstrainedMap const& aEdge_is_constrained_map // defaults to No_constrained_edge_map<TM>()
const VertexIndexMap& aVertex_index_map, // defaults to get(vertex_index, aSurface)
const VertexPointMap& aVertex_point_map, // defaults to get(vertex_point, aSurface)
const EdgeIndexMap& aEdge_index_map, // defaults to get(edge_index, aSurface)
const EdgeIsConstrainedMap& aEdge_is_constrained_map, // defaults to No_constrained_edge_map<TM>()
// optional strategy policies - defaults to LindstomTurk
, GetCost const& aGet_cost
, GetPlacement const& aGet_placement
, Visitor aVisitor
)
const GetCost& aGet_cost,
const GetPlacement& aGet_placement,
Visitor aVisitor)
{
typedef EdgeCollapse< TM
, ShouldStop
, VertexIndexMap
, VertexPointMap
, EdgeIndexMap
, EdgeIsConstrainedMap
, GetCost
, GetPlacement
, Visitor
>
Algorithm;
typedef EdgeCollapse<TM, ShouldStop,
VertexIndexMap, VertexPointMap, EdgeIndexMap, EdgeIsConstrainedMap,
GetCost, GetPlacement, Visitor> Algorithm;
Algorithm algorithm( aSurface
, aShould_stop
, aVertex_index_map
, aVertex_point_map
, aEdge_index_map
, aEdge_is_constrained_map
, aGet_cost
, aGet_placement
, aVisitor
) ;
Algorithm algorithm(aSurface, aShould_stop,
aVertex_index_map, aVertex_point_map, aEdge_index_map, aEdge_is_constrained_map,
aGet_cost, aGet_placement, aVisitor);
return algorithm.run();
}
struct Dummy_visitor
{
template<class TM> void OnStarted( TM& ) const {}
template<class TM> void OnFinished ( TM& ) const {}
template<class Profile> void OnStopConditionReached( Profile const& ) const {}
template<class Profile, class OFT> void OnCollected( Profile const&, OFT const& ) const {}
template<class Profile, class OFT, class Size_type> void OnSelected( Profile const&, OFT const&, Size_type, Size_type ) const {}
template<class Profile, class OPoint> void OnCollapsing(Profile const&, OPoint const& ) const {}
template<class Profile, class VH> void OnCollapsed( Profile const&, VH ) const {}
template<class Profile> void OnNonCollapsable(Profile const& ) const {}
} ;
template<class TM>
void OnStarted(TM&) const {}
template<class TM>
void OnFinished(TM&) const {}
template<class Profile>
void OnStopConditionReached(const Profile&) const {}
template<class Profile, class OFT>
void OnCollected(Profile const&, const OFT&) const {}
template<class Profile, class OFT, class Size_type>
void OnSelected(Profile const&, const OFT&, Size_type, Size_type) const {}
template<class Profile, class OPoint>
void OnCollapsing(const Profile&, const OPoint&) const {}
template<class Profile, class VH>
void OnCollapsed(const Profile&, VH) const {}
template<class Profile>
void OnNonCollapsable(Profile const&) const {}
};
template<class TM, class ShouldStop, class P, class T, class R>
int edge_collapse ( TM& aSurface
, ShouldStop const& aShould_stop
, cgal_bgl_named_params<P,T,R> const& aParams
)
int edge_collapse(TM& aSurface,
const ShouldStop& aShould_stop,
cgal_bgl_named_params<P,T,R> const& aParams)
{
using boost::choose_param ;
using boost::choose_const_pmap ;
using boost::get_param ;
using boost::choose_param;
using boost::choose_const_pmap;
using boost::get_param;
LindstromTurk_params lPolicyParams ;
internal_np::graph_visitor_t vis = internal_np::graph_visitor_t() ;
return edge_collapse(aSurface
,aShould_stop
,choose_const_pmap(get_param(aParams,internal_np::vertex_index),aSurface,boost::vertex_index)
,choose_pmap(get_param(aParams,internal_np::vertex_point),aSurface,boost::vertex_point)
,choose_const_pmap(get_param(aParams,internal_np::halfedge_index),aSurface,boost::halfedge_index)
,choose_param (get_param(aParams,internal_np::edge_is_constrained),No_constrained_edge_map<TM>())
,choose_param (get_param(aParams,internal_np::get_cost_policy), LindstromTurk_cost<TM>())
,choose_param (get_param(aParams,internal_np::get_placement_policy), LindstromTurk_placement<TM>())
,choose_param (get_param(aParams,vis), Dummy_visitor())
);
LindstromTurk_params lPolicyParams;
internal_np::graph_visitor_t vis = internal_np::graph_visitor_t();
return edge_collapse(aSurface, aShould_stop,
choose_const_pmap(get_param(aParams,internal_np::vertex_index),aSurface,boost::vertex_index),
choose_pmap(get_param(aParams,internal_np::vertex_point),aSurface,boost::vertex_point),
choose_const_pmap(get_param(aParams,internal_np::halfedge_index),aSurface,boost::halfedge_index),
choose_param(get_param(aParams,internal_np::edge_is_constrained),No_constrained_edge_map<TM>()),
choose_param(get_param(aParams,internal_np::get_cost_policy), LindstromTurk_cost<TM>()),
choose_param(get_param(aParams,internal_np::get_placement_policy), LindstromTurk_placement<TM>()),
choose_param(get_param(aParams,vis), Dummy_visitor()));
}
template<class TM, class ShouldStop, class GT, class P, class T, class R>
int edge_collapse ( TM& aSurface
, ShouldStop const& aShould_stop
, cgal_bgl_named_params<P,T,R> const& aParams
)
template<class TM, class ShouldStop, class GT, class P, class T, class R>
int edge_collapse(TM& aSurface,
ShouldStop const& aShould_stop,
cgal_bgl_named_params<P,T,R> const& aParams)
{
using boost::choose_param ;
using boost::choose_const_pmap ;
using boost::get_param ;
using boost::choose_param;
using boost::choose_const_pmap;
using boost::get_param;
LindstromTurk_params lPolicyParams ;
internal_np::graph_visitor_t vis = internal_np::graph_visitor_t() ;
return edge_collapse(aSurface
,aShould_stop
,choose_const_pmap(get_param(aParams,internal_np::vertex_index),aSurface,boost::vertex_index)
,choose_const_pmap(get_param(aParams,internal_np::vertex_point),aSurface,boost::vertex_point)
,choose_const_pmap(get_param(aParams,internal_np::halfedge_index),aSurface,boost::halfedge_index)
,choose_param (get_param(aParams,internal_np::edge_is_constrained),No_constrained_edge_map<TM>())
,choose_param (get_param(aParams,internal_np::get_cost_policy), LindstromTurk_cost<TM>())
,choose_param (get_param(aParams,internal_np::get_placement_policy), LindstromTurk_placement<TM>())
,choose_param (get_param(aParams,vis), Dummy_visitor())
);
LindstromTurk_params lPolicyParams;
internal_np::graph_visitor_t vis = internal_np::graph_visitor_t();
return edge_collapse(aSurface, aShould_stop,
choose_const_pmap(get_param(aParams,internal_np::vertex_index),aSurface,boost::vertex_index),
choose_const_pmap(get_param(aParams,internal_np::vertex_point),aSurface,boost::vertex_point),
choose_const_pmap(get_param(aParams,internal_np::halfedge_index),aSurface,boost::halfedge_index),
choose_param(get_param(aParams,internal_np::edge_is_constrained),No_constrained_edge_map<TM>()),
choose_param(get_param(aParams,internal_np::get_cost_policy), LindstromTurk_cost<TM>()),
choose_param(get_param(aParams,internal_np::get_placement_policy), LindstromTurk_placement<TM>()),
choose_param(get_param(aParams,vis), Dummy_visitor()));
}
template<class TM, class ShouldStop>
int edge_collapse ( TM& aSurface, ShouldStop const& aShould_stop )
int edge_collapse(TM& aSurface, const ShouldStop& aShould_stop)
{
return edge_collapse(aSurface,aShould_stop, CGAL::parameters::halfedge_index_map(get(boost::halfedge_index,aSurface)));
}
template<class TM, class ShouldStop, class GT>
int edge_collapse ( TM& aSurface, ShouldStop const& aShould_stop)
template<class TM, class ShouldStop, class GT>
int edge_collapse(TM& aSurface, const ShouldStop& aShould_stop)
{
return edge_collapse(aSurface,aShould_stop, CGAL::parameters::halfedge_index_map(get(boost::halfedge_index,aSurface)));
}
} // namespace Surface_mesh_simplification
} //namespace CGAL
#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_EDGE_COLLAPSE_H //
// EOF //
} // namespace CGAL
#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_EDGE_COLLAPSE_H

16
Surface_mesh_simplification/test/Surface_mesh_simplification/basics.h
View File

@ -16,10 +16,10 @@
//#define CGAL_SURFACE_SIMPLIFICATION_ENABLE_TRACE 4
//#define CGAL_SURFACE_SIMPLIFICATION_ENABLE_LT_TRACE 4
void Surface_simplification_external_trace( std::string s )
void Surface_simplification_external_trace(std::string s)
{
static std::ofstream out("log.txt");
out << s << std::endl ;
out << s << std::endl;
}
#include <CGAL/Real_timer.h>
@ -39,9 +39,9 @@ void Surface_simplification_external_trace( std::string s )
#include <CGAL/assertions_behaviour.h>
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel ;
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::FT FT ;
typedef Kernel::FT FT;
typedef Kernel::Vector_3 Vector;
typedef Kernel::Point_3 Point;
@ -60,11 +60,11 @@ typedef Surface::Facet_const_iterator Facet_const_iterator;
typedef Surface::Facet_const_handle Facet_const_handle;
typedef Surface::Halfedge_around_vertex_const_circulator HV_circulator;
typedef Surface::Halfedge_around_facet_circulator HF_circulator;
typedef Surface::size_type size_type ;
typedef Surface::size_type size_type;
using namespace std ;
using namespace boost ;
using namespace CGAL ;
using namespace std;
using namespace boost;
using namespace CGAL;

20
Surface_mesh_simplification/test/Surface_mesh_simplification/edge_collapse_topology.cpp
View File

@ -13,27 +13,27 @@
typedef CGAL::Simple_cartesian<double> Kernel;
typedef CGAL::Polyhedron_3<Kernel> Surface;
namespace SMS = CGAL::Surface_mesh_simplification ;
namespace SMS = CGAL::Surface_mesh_simplification;
int main( int argc, char** argv )
int main(int argc, char** argv)
{
if (argc!=2){
if(argc!=2){
std::cerr << "Please provide only an off-file as input\n";
return 1;
}
std::ifstream is(argv[1]);
if (!is){
if(!is){
std::cerr << "Error reading the input\n";
return 1;
}
Surface surface;
is >> surface ;
is >> surface;
std::size_t initial_count = (surface.size_of_halfedges()/2);
std::cout << "Initial count " << initial_count << " edges.\n" ;
std::cout << "Initial count " << initial_count << " edges.\n";
// Contract the surface as much as possible
SMS::Count_stop_predicate<Surface> stop(0);
@ -46,12 +46,12 @@ int main( int argc, char** argv )
);
std::cout << "\nFinished...\n" << r << " edges removed.\n"
<< (surface.size_of_halfedges()/2) << " final edges.\n" ;
<< (surface.size_of_halfedges()/2) << " final edges.\n";
assert( initial_count == (surface.size_of_halfedges()/2) + r );
// std::ofstream os( argc > 2 ? argv[2] : "out.off" ) ; os << surface ;
assert(initial_count == (surface.size_of_halfedges()/2) + r);
// std::ofstream os(argc > 2 ? argv[2] : "out.off"); os << surface;
return 0 ;
return 0;
}
// EOF //

328
Surface_mesh_simplification/test/Surface_mesh_simplification/test_edge_collapse_Polyhedron_3.cpp
View File

@ -28,10 +28,10 @@
//#define TEST_TEST_TRACE_ENABLED
template<std::size_t N> struct eat_sizeof {} ;
template<std::size_t N> struct eat_sizeof {};
#ifdef TEST_TEST_TRACE_ENABLED
# define TEST_TRACE(m) std::cerr << m << std::endl ;
# define TEST_TRACE(m) std::cerr << m << std::endl;
#else
# define TEST_TRACE(m) (eat_sizeof<sizeof(m)>())
@ -39,54 +39,54 @@ template<std::size_t N> struct eat_sizeof {} ;
// This is here only to allow a breakpoint to be placed so I can trace back the problem.
void error_handler ( char const* what, char const* expr, char const* file, int line, char const* msg )
void error_handler (char const* what, char const* expr, char const* file, int line, char const* msg)
{
cerr << "CGAL error: " << what << " violation!" << endl
<< "Expr: " << expr << endl
<< "File: " << file << endl
<< "Line: " << line << endl;
if ( msg != 0)
if(msg != 0)
cerr << "Explanation:" << msg << endl;
throw std::runtime_error(expr);
}
namespace SMS = CGAL::Surface_mesh_simplification ;
namespace SMS = CGAL::Surface_mesh_simplification;
typedef SMS::Edge_profile<Surface> Profile ;
typedef SMS::Edge_profile<Surface> Profile;
typedef boost::shared_ptr<Surface> SurfaceSP ;
typedef boost::shared_ptr<Surface> SurfaceSP;
// Constructs a flat polyhedron containing just the link of an edge or vertex.
class Link_builder : public CGAL::Modifier_base<Surface::HalfedgeDS>
{
map<std::size_t, std::size_t> mMap ;
map<std::size_t, std::size_t> mMap;
int mVIdx ;
int mVIdx;
protected:
typedef boost::graph_traits<Surface> GraphTraits ;
typedef boost::graph_traits<Surface> GraphTraits;
typedef CGAL::Halfedge_around_source_iterator<Surface> out_edge_iterator ;
typedef GraphTraits::halfedge_descriptor halfedge_descriptor ;
typedef GraphTraits::vertex_descriptor vertex_descriptor ;
typedef CGAL::Polyhedron_incremental_builder_3<Surface::HalfedgeDS> Builder ;
typedef CGAL::Halfedge_around_source_iterator<Surface> out_edge_iterator;
typedef GraphTraits::halfedge_descriptor halfedge_descriptor;
typedef GraphTraits::vertex_descriptor vertex_descriptor;
typedef CGAL::Polyhedron_incremental_builder_3<Surface::HalfedgeDS> Builder;
Link_builder() : mVIdx(0) {}
void add_vertex( Builder& aBuilder, Vertex_const_handle v )
void add_vertex(Builder& aBuilder, Vertex_const_handle v)
{
aBuilder.add_vertex(v->point());
mMap.insert( make_pair(v->id(),mVIdx++) );
mMap.insert(make_pair(v->id(),mVIdx++));
}
void add_triangle( Builder& aBuilder, Profile::Triangle const& aTri )
void add_triangle(Builder& aBuilder, Profile::Triangle const& aTri)
{
aBuilder.begin_facet();
aBuilder.add_vertex_to_facet( mMap[aTri.v0->id()] );
aBuilder.add_vertex_to_facet( mMap[aTri.v1->id()] );
aBuilder.add_vertex_to_facet( mMap[aTri.v2->id()] );
aBuilder.add_vertex_to_facet(mMap[aTri.v0->id()]);
aBuilder.add_vertex_to_facet(mMap[aTri.v1->id()]);
aBuilder.add_vertex_to_facet(mMap[aTri.v2->id()]);
aBuilder.end_facet();
}
};
@ -94,26 +94,26 @@ protected:
// Constructs a flat polyhedron containing just the link of an edge
class Edge_link_builder : public Link_builder
{
Profile const& mProfile ;
const Profile& mProfile;
public:
Edge_link_builder( Profile const& aProfile ) : mProfile(aProfile) {}
Edge_link_builder(const Profile& aProfile) : mProfile(aProfile) {}
void operator()( Surface::HalfedgeDS& hds)
void operator()(Surface::HalfedgeDS& hds)
{
Builder B( hds, true);
Builder B(hds, true);
B.begin_surface( 2 + mProfile.link().size(), mProfile.triangles().size() );
B.begin_surface(2 + mProfile.link().size(), mProfile.triangles().size());
this->add_vertex(B,mProfile.v0());
this->add_vertex(B,mProfile.v1());
for ( Profile::vertex_descriptor_vector::const_iterator vit = mProfile.link().begin(); vit != mProfile.link().end(); ++ vit )
this->add_vertex(B, *vit );
for(Profile::vertex_descriptor_vector::const_iterator vit = mProfile.link().begin(); vit != mProfile.link().end(); ++ vit)
this->add_vertex(B, *vit);
for ( Profile::Triangle_vector::const_iterator tit = mProfile.triangles().begin(); tit != mProfile.triangles().end(); ++ tit )
for(Profile::Triangle_vector::const_iterator tit = mProfile.triangles().begin(); tit != mProfile.triangles().end(); ++ tit)
this->add_triangle(B,*tit);
B.end_surface();
@ -123,227 +123,227 @@ public:
// Constructs a flat polyhedron containing just the link of a vertex
class Vertex_link_builder : public Link_builder
{
Surface* m_tm ;
Vertex_handle mV ;
Surface* m_tm;
Vertex_handle mV;
Surface& tm() { return *m_tm ; }
Surface& tm() { return *m_tm; }
public:
Vertex_link_builder( Surface& aTM, Vertex_handle aV) : m_tm(&aTM), mV(aV) {}
Vertex_link_builder(Surface& aTM, Vertex_handle aV) : m_tm(&aTM), mV(aV) {}
void operator()( Surface::HalfedgeDS& hds )
void operator()(Surface::HalfedgeDS& hds)
{
Builder B( hds, true);
Builder B(hds, true);
B.begin_surface( 1 + out_degree(mV,tm()), out_degree(mV,tm()) );
B.begin_surface(1 + out_degree(mV,tm()), out_degree(mV,tm()));
this->add_vertex(B,mV);
Profile::Triangle_vector triangles ;
Profile::Triangle_vector triangles;
out_edge_iterator eb, ee ;
for ( boost::tie(eb,ee) = halfedges_around_source(opposite(halfedge(mV,tm()),tm()),tm()) ; eb != ee ; ++ eb )
out_edge_iterator eb, ee;
for(boost::tie(eb,ee) = halfedges_around_source(opposite(halfedge(mV,tm()),tm()),tm()); eb != ee; ++ eb)
{
halfedge_descriptor out_edge1 = *eb ;
halfedge_descriptor out_edge1 = *eb;
halfedge_descriptor out_edge2 = out_edge1->opposite()->next();
vertex_descriptor v1 = target(out_edge1,tm());
vertex_descriptor v2 = target(out_edge2,tm());
this->add_vertex(B,v1);
triangles.push_back( Profile::Triangle(mV,v1,v2) ) ;
triangles.push_back(Profile::Triangle(mV,v1,v2));
}
for ( Profile::Triangle_vector::const_iterator tit = triangles.begin(); tit != triangles.end(); ++ tit )
for(Profile::Triangle_vector::const_iterator tit = triangles.begin(); tit != triangles.end(); ++ tit)
this->add_triangle(B,*tit);
B.end_surface();
}
};
SurfaceSP create_edge_link ( Profile const& aProfile )
SurfaceSP create_edge_link (const Profile& aProfile)
{
SurfaceSP rLink( new Surface );
SurfaceSP rLink(new Surface);
try
{
Edge_link_builder lBuilder(aProfile) ;
Edge_link_builder lBuilder(aProfile);
rLink->delegate(lBuilder);
}
catch(...)
{
cerr << "Unable to create local mesh" << endl ;
cerr << "Unable to create local mesh" << endl;
}
return rLink ;
return rLink;
}
SurfaceSP create_vertex_link ( Profile const& aProfile, Vertex_handle aV )
SurfaceSP create_vertex_link (const Profile& aProfile, Vertex_handle aV)
{
SurfaceSP rLink( new Surface );
SurfaceSP rLink(new Surface);
try
{
Vertex_link_builder lBuilder(aProfile.surface(),aV) ;
Vertex_link_builder lBuilder(aProfile.surface(),aV);
rLink->delegate(lBuilder);
}
catch ( ... )
catch (...)
{
cerr << "Unable to create local mesh" << endl ;
cerr << "Unable to create local mesh" << endl;
}
return rLink ;
return rLink;
}
void write ( SurfaceSP aSurface, string aName )
void write (SurfaceSP aSurface, string aName)
{
ofstream out(aName.c_str());
out << *aSurface ;
out << *aSurface;
}
template<class T>
string opt2str ( boost::optional<T> const& o )
string opt2str (boost::optional<T> const& o)
{
ostringstream ss ;
if ( o )
ss << *o ;
else ss << "<none>" ;
ostringstream ss;
if(o)
ss << *o;
else ss << "<none>";
return ss.str();
}
template<class N>
string float2str ( N n )
string float2str (N n)
{
return boost::str( boost::format("%+05.19g") % to_double(n) ) ;
return boost::str(boost::format("%+05.19g") % to_double(n));
}
template<class P>
string point2str ( P const& p )
string point2str (P const& p)
{
return boost::str( boost::format( "(%+05.19g %+05.19g %+05.19g)") % to_double(p.x()) % to_double(p.y()) % to_double(p.z()) ) ;
return boost::str(boost::format("(%+05.19g %+05.19g %+05.19g)") % to_double(p.x()) % to_double(p.y()) % to_double(p.z()));
}
template<class P>
string optpoint2str ( boost::optional<P> const& p )
string optpoint2str (boost::optional<P> const& p)
{
ostringstream ss ;
if ( p )
ostringstream ss;
if(p)
ss << point2str(*p);
else ss << "<none>" ;
else ss << "<none>";
return ss.str();
}
template<class N>
string optfloat2str ( boost::optional<N> const& n )
string optfloat2str (boost::optional<N> const& n)
{
ostringstream ss ;
if ( n )
ostringstream ss;
if(n)
ss << float2str(*n);
else ss << "<none>" ;
else ss << "<none>";
return ss.str();
}
template<class V>
string vertex2str ( V const& v )
string vertex2str (V const& v)
{
ostringstream ss ;
ss << "[V" << v->id() << point2str(v->point()) << "]" ;
ostringstream ss;
ss << "[V" << v->id() << point2str(v->point()) << "]";
return ss.str();
}
template<class E>
string edge2str ( E const& e )
string edge2str (E const& e)
{
ostringstream ss ;
ss << "{E" << e->id() << vertex2str(e->opposite()->vertex()) << "->" << vertex2str(e->vertex()) << "}" ;
ostringstream ss;
ss << "{E" << e->id() << vertex2str(e->opposite()->vertex()) << "->" << vertex2str(e->vertex()) << "}";
return ss.str();
}
template<class T> ostream& operator << ( ostream& os, boost::optional<T> const& o ) { return os << opt2str(o); }
template<class T> ostream& operator << (ostream& os, boost::optional<T> const& o) { return os << opt2str(o); }
string normalize_EOL ( string line )
string normalize_EOL (string line)
{
string::size_type l = line.length();
string::size_type d = ( l > 0 && line[l-1] == '\r' ) ? 1 : 0 ;
return line.substr(0, l-d ) ;
string::size_type d = (l > 0 && line[l-1] == '\r') ? 1 : 0;
return line.substr(0, l-d);
}
#define REPORT_ERROR(msg) error(__FILE__,__LINE__,0, msg);
#define REPORT_ERROR2(pred,msg) REPORT_ERROR(msg)
#define CHECK_MSG(pred,msg) if (!(pred)) REPORT_ERROR2(#pred,msg)
#define CHECK_MSG(pred,msg) if(!(pred)) REPORT_ERROR2(#pred,msg)
#define CHECK(pred) CHECK_MSG(pred,string(""))
#define CHECK_EQUAL(x,y) CHECK_MSG(((x)==(y)), str(format("Assertion failed: %1%(=%2%)==%3%(=%4%)") % (#x) % (x) % (#y) % (y) ) )
#define CHECK_NOT_EQUAL(x,y) CHECK_MSG(((x)!=(y)), str(format("Assertion failed: %1%(=%2%)!=%3%(=%4%)") % (#x) % (x) % (#y) % (y) ) )
#define CHECK_EQUAL(x,y) CHECK_MSG(((x)==(y)), str(format("Assertion failed: %1%(=%2%)==%3%(=%4%)") % (#x) % (x) % (#y) % (y)))
#define CHECK_NOT_EQUAL(x,y) CHECK_MSG(((x)!=(y)), str(format("Assertion failed: %1%(=%2%)!=%3%(=%4%)") % (#x) % (x) % (#y) % (y)))
class Visitor : public SMS::Edge_collapse_visitor_base<Surface>
{
public :
Visitor( string aTestCase ) : mTestCase(aTestCase)
Visitor(string aTestCase) : mTestCase(aTestCase)
{
#ifdef CGAL_SMS_TRACE_IMPL
::internal::cgal_enable_sms_trace = true ;
::internal::cgal_enable_sms_trace = true;
#endif
mStep = 0 ;
mStep = 0;
}
virtual ~Visitor()
{
}
virtual void OnStarted( Surface& ) const
virtual void OnStarted(Surface&) const
{
}
virtual void OnFinished ( Surface& aSurface ) const
virtual void OnFinished (Surface& aSurface) const
{
CHECK(aSurface.is_valid());
}
virtual void OnCollected( Profile const& aProfile, boost::optional<FT> const& aCost ) const
virtual void OnCollected(const Profile& aProfile, boost::optional<FT> const& aCost) const
{
TEST_TRACE( str ( format("Collecting %1% : cost=%2%") % edge2str(aProfile.v0_v1()) % optfloat2str(aCost) ) ) ;
TEST_TRACE(str (format("Collecting %1% : cost=%2%") % edge2str(aProfile.v0_v1()) % optfloat2str(aCost)));
}
virtual void OnCollapsing( Profile const& aProfile, boost::optional<Point> const& aP ) const
virtual void OnCollapsing(const Profile& aProfile, boost::optional<Point> const& aP) const
{
TEST_TRACE( str ( format("S %1% - Collapsing %2% : placement=%3%") % mStep % edge2str(aProfile.v0_v1()) % optpoint2str(aP) ) ) ;
TEST_TRACE(str (format("S %1% - Collapsing %2% : placement=%3%") % mStep % edge2str(aProfile.v0_v1()) % optpoint2str(aP)));
//mBefore = create_edge_link(aProfile);
}
virtual void OnCollapsed( Profile const& aProfile, Vertex_handle const& aV ) const
virtual void OnCollapsed(const Profile& aProfile, Vertex_handle const& aV) const
{
// Some collapse can result in self-intersections in any case, so we can't mark it as a failure
if ( false && Is_self_intersecting( aProfile.surface() ) )
if(false && Is_self_intersecting(aProfile.surface()))
{
SurfaceSP lAfter = create_vertex_link(aProfile, aV);
write(mBefore, str( format( "%1%.step-%2%-before.off" ) % mTestCase % mStep ) ) ;
write(lAfter , str( format( "%1%.step-%2%-after.off" ) % mTestCase % mStep ) ) ;
write(mBefore, str(format("%1%.step-%2%-before.off") % mTestCase % mStep));
write(lAfter , str(format("%1%.step-%2%-after.off") % mTestCase % mStep));
REPORT_ERROR( str( format("Resulting surface self-intersects after step %1% (%2% edges left)") % mStep % ( aProfile.surface().size_of_halfedges() / 2 ) ) ) ;
REPORT_ERROR(str(format("Resulting surface self-intersects after step %1% (%2% edges left)") % mStep % (aProfile.surface().size_of_halfedges() / 2)));
}
++ mStep ;
++ mStep;
}
private:
void error ( char const* file, int line, char const* pred, string msg ) const
void error (char const* file, int line, char const* pred, string msg) const
{
cerr << "ERROR in " << file << " at " << line << endl ;
if ( pred )
cerr << " Assertion failed: " << pred << endl ;
cerr << " " << msg << endl ;
cerr << "ERROR in " << file << " at " << line << endl;
if(pred)
cerr << " Assertion failed: " << pred << endl;
cerr << " " << msg << endl;
throw runtime_error("");
@ -351,100 +351,100 @@ private:
private:
string mTestCase ;
mutable int mStep ;
string mTestCase;
mutable int mStep;
mutable SurfaceSP mBefore;
} ;
};
bool sWriteResult = false ;
bool sUseMP = false ;
int sStop = 1 ;
bool sWriteResult = false;
bool sUseMP = false;
int sStop = 1;
bool Test ( string aName )
bool Test (string aName)
{
bool rSucceeded = false ;
bool rSucceeded = false;
try
{
string off_name = aName ;
string off_name = aName;
ifstream off_is(off_name.c_str());
if ( off_is )
if(off_is)
{
Surface lSurface;
off_is >> lSurface ;
if ( lSurface.is_valid() )
off_is >> lSurface;
if(lSurface.is_valid())
{
if ( lSurface.is_pure_triangle() )
if(lSurface.is_pure_triangle())
{
cerr << "Processing " << aName << " (" << ( lSurface.size_of_halfedges() / 2 ) << " edges)" << endl ;
cerr << "Processing " << aName << " (" << (lSurface.size_of_halfedges() / 2) << " edges)" << endl;
Visitor vis(aName) ;
Visitor vis(aName);
set_halfedgeds_items_id(lSurface);
SMS::Count_stop_predicate<Surface> stop(sStop);
Real_timer t ; t.start();
Real_timer t; t.start();
if ( sUseMP )
if(sUseMP)
{
SMS::Edge_length_cost <Surface> cost ;
SMS::Midpoint_placement<Surface> placement ;
SMS::Edge_length_cost <Surface> cost;
SMS::Midpoint_placement<Surface> placement;
edge_collapse(lSurface,stop,CGAL::parameters::get_cost(cost).get_placement(placement).visitor(vis) );
edge_collapse(lSurface,stop,CGAL::parameters::get_cost(cost).get_placement(placement).visitor(vis));
}
else
{
SMS::LindstromTurk_cost <Surface> cost ;
SMS::LindstromTurk_placement<Surface> placement ;
SMS::LindstromTurk_cost <Surface> cost;
SMS::LindstromTurk_placement<Surface> placement;
edge_collapse(lSurface,stop,CGAL::parameters::get_cost(cost).get_placement(placement).visitor(vis) );
edge_collapse(lSurface,stop,CGAL::parameters::get_cost(cost).get_placement(placement).visitor(vis));
}
t.stop();
rSucceeded = lSurface.is_valid() && lSurface.is_pure_triangle() ;
rSucceeded = lSurface.is_valid() && lSurface.is_pure_triangle();
cerr << "\r" << aName << ( rSucceeded ? " succeeded" : "FAILED" ) << ". Elapsed time=" << t.time() << " seconds." << endl ;
cerr << "\r" << aName << (rSucceeded ? " succeeded" : "FAILED") << ". Elapsed time=" << t.time() << " seconds." << endl;
if ( sWriteResult )
if(sWriteResult)
{
string out_name = off_name ;
out_name.replace( off_name.find_last_of('.'), 15, ".simplified.off");
string out_name = off_name;
out_name.replace(off_name.find_last_of('.'), 15, ".simplified.off");
ofstream out(out_name.c_str());
out << lSurface ;
cerr << "Result written into: " << out_name << endl ;
out << lSurface;
cerr << "Result written into: " << out_name << endl;
}
}
else
{
cerr << "Surface is not triangulated (has faces with more than 3 sides): " << aName << endl ;
cerr << "Surface is not triangulated (has faces with more than 3 sides): " << aName << endl;
}
}
else
{
cerr << "Invalid surface: " << aName << endl ;
cerr << "Invalid surface: " << aName << endl;
}
}
else
{
cerr << "Unable to open test file " << aName << endl ;
cerr << "Unable to open test file " << aName << endl;
}
}
catch ( std::exception const& x )
catch (std::exception const& x)
{
string what(x.what());
if ( what.length() > 0 )
cerr << "Exception caught: " << what << endl ;
if(what.length() > 0)
cerr << "Exception caught: " << what << endl;
}
return rSucceeded ;
return rSucceeded;
}
int main( int argc, char** argv )
int main(int argc, char** argv)
{
set_error_handler (error_handler);
set_warning_handler(error_handler);
@ -452,52 +452,52 @@ int main( int argc, char** argv )
cout << setprecision(4);
cerr << setprecision(4);
vector<string> lCases ;
vector<string> lCases;
for ( int i = 1 ; i < argc ; ++i )
for(int i = 1; i < argc; ++i)
{
if ( argv[i][0] == '-' )
if(argv[i][0] == '-')
{
switch(argv[i][1])
{
case 'w': sWriteResult = true ; break ;
case 's': sStop = atoi(&argv[i][2]) ; break ;
case 'm': sUseMP = true ; break ;
case 'w': sWriteResult = true; break;
case 's': sStop = atoi(&argv[i][2]); break;
case 'm': sUseMP = true; break;
}
}
else
{
string c( normalize_EOL( string(argv[i]) ) ) ;
string::size_type pos = c.find_last_of(".") ;
if ( pos != string::npos )
string c(normalize_EOL(string(argv[i])));
string::size_type pos = c.find_last_of(".");
if(pos != string::npos)
{
string ext = c.substr(pos);
if ( ext == ".off" )
if(ext == ".off")
lCases.push_back(c);
}
}
}
if ( lCases.size() == 0 )
if(lCases.size() == 0)
{
cout << "collapse_edge_test file0 [-w] [-sNUM_EDGES] [-m] file1 ... fileN\n" ;
return 1 ;
cout << "collapse_edge_test file0 [-w] [-sNUM_EDGES] [-m] file1 ... fileN\n";
return 1;
}
else
{
unsigned lOK = 0 ;
for ( vector<string>::const_iterator it = lCases.begin(); it != lCases.end() ; ++ it )
unsigned lOK = 0;
for(vector<string>::const_iterator it = lCases.begin(); it != lCases.end(); ++ it)
{
if ( Test(*it) )
++ lOK ;
if(Test(*it))
++ lOK;
}
cout << endl
<< lOK << " cases succedded." << endl
<< (lCases.size() - lOK ) << " cases failed." << endl ;
<< (lCases.size() - lOK) << " cases failed." << endl;
return lOK == lCases.size() ? 0 : 1 ;
return lOK == lCases.size() ? 0 : 1;
}
}
@ -505,15 +505,15 @@ int main( int argc, char** argv )
namespace CGAL
{
void assertion_fail( const char* expr, const char* file, int line )
void assertion_fail(const char* expr, const char* file, int line)
{
assertion_fail(expr,file,line,"");
}
void precondition_fail( const char* expr, const char* file, int line )
void precondition_fail(const char* expr, const char* file, int line)
{
precondition_fail(expr,file,line,"");
}
void postcondition_fail( const char* expr, const char* file, int line )
void postcondition_fail(const char* expr, const char* file, int line)
{
postcondition_fail(expr,file,line,"");
}

28
Surface_mesh_simplification/test/Surface_mesh_simplification/test_edge_profile_link.cpp
View File

@ -13,12 +13,12 @@ void naive_all_triangles(Mesh::Halfedge_index h, Mesh& m, std::set<Mesh::Face_in
{
for(Mesh::Halfedge_index hh : CGAL::halfedges_around_source(h, m))
{
if (!is_border(hh, m))
if(!is_border(hh, m))
triangles.insert(face(hh,m));
}
for(Mesh::Halfedge_index hh : CGAL::halfedges_around_target(h, m))
{
if (!is_border(hh, m))
if(!is_border(hh, m))
triangles.insert(face(hh,m));
}
}
@ -34,8 +34,8 @@ void naive_link_vertices(Mesh::Halfedge_index h, Mesh& m,
link_vertices.insert(target(h, m));
}
}
link_vertices.erase( source(h, m) );
link_vertices.erase( target(h, m) );
link_vertices.erase(source(h, m));
link_vertices.erase(target(h, m));
}
struct A{};
@ -65,7 +65,7 @@ void test(const char* fname)
{
Mesh m;
std::ifstream input(fname);
assert( !input.fail() );
assert(!input.fail());
input >> m;
assert(num_vertices(m)!=0);
A a;
@ -77,18 +77,18 @@ void test(const char* fname)
Profile profile(h, m, a, get(boost::vertex_point, m), a, true);
if (CGAL::Euler::does_satisfy_link_condition(edge(h, m), m))
if(CGAL::Euler::does_satisfy_link_condition(edge(h, m), m))
{
std::set<Mesh::Vertex_index> link_vertices;
naive_link_vertices(h, m, triangles, link_vertices);
assert( link_vertices.size()==profile.link().size() );
assert( std::set<Mesh::Vertex_index>(profile.link().begin(),
assert(link_vertices.size()==profile.link().size());
assert(std::set<Mesh::Vertex_index>(profile.link().begin(),
profile.link().end()).size()
== link_vertices.size() );
== link_vertices.size());
for(const Mesh::Vertex_index& v : profile.link())
{
assert( link_vertices.count(v) == 1 );
assert(link_vertices.count(v) == 1);
}
}
@ -96,23 +96,23 @@ void test(const char* fname)
std::set<boost::tuple<Mesh::Vertex_index, Mesh::Vertex_index, Mesh::Vertex_index> > triple_set;
for(const Profile::Triangle& t : profile.triangles())
{
triple_set.insert( make_canonical_tuple(t) );
triple_set.insert(make_canonical_tuple(t));
}
for(Mesh::Face_index f : triangles)
{
assert( triple_set.count( make_canonical_tuple(f, m) ) == 1 );
assert(triple_set.count(make_canonical_tuple(f, m)) == 1);
}
}
}
int main(int argc, char** argv)
{
for (int i=1; i<argc; ++i)
for(int i=1; i<argc; ++i)
{
std::cout << "Testing " << argv[i] << "\n";
test(argv[i]);
}
if (argc==1)
if(argc==1)
{
std::cout << "No file provided, nothing tested\n";
}

66
Surface_mesh_simplification/test/Surface_mesh_simplification/test_self_intersection.h
View File

@ -13,61 +13,61 @@ std::vector<Triangle> triangles;
struct Intersect_facets
{
void operator()( const Box* b, const Box* c) const
void operator()(const Box* b, const Box* c) const
{
Halfedge_const_handle h = b->handle()->halfedge();
// check for shared egde --> no intersection
if ( h->opposite()->facet() == c->handle()
if(h->opposite()->facet() == c->handle()
|| h->next()->opposite()->facet() == c->handle()
|| h->next()->next()->opposite()->facet() == c->handle())
return;
// check for shared vertex --> maybe intersection, maybe not
Halfedge_const_handle g = c->handle()->halfedge();
Halfedge_const_handle v;
if ( h->vertex() == g->vertex())
if(h->vertex() == g->vertex())
v = g;
if ( h->vertex() == g->next()->vertex())
if(h->vertex() == g->next()->vertex())
v = g->next();
if ( h->vertex() == g->next()->next()->vertex())
if(h->vertex() == g->next()->next()->vertex())
v = g->next()->next();
if ( v == Halfedge_const_handle()) {
if(v == Halfedge_const_handle()) {
h = h->next();
if ( h->vertex() == g->vertex())
if(h->vertex() == g->vertex())
v = g;
if ( h->vertex() == g->next()->vertex())
if(h->vertex() == g->next()->vertex())
v = g->next();
if ( h->vertex() == g->next()->next()->vertex())
if(h->vertex() == g->next()->next()->vertex())
v = g->next()->next();
if ( v == Halfedge_const_handle()) {
if(v == Halfedge_const_handle()) {
h = h->next();
if ( h->vertex() == g->vertex())
if(h->vertex() == g->vertex())
v = g;
if ( h->vertex() == g->next()->vertex())
if(h->vertex() == g->next()->vertex())
v = g->next();
if ( h->vertex() == g->next()->next()->vertex())
if(h->vertex() == g->next()->next()->vertex())
v = g->next()->next();
}
}
if ( v != Halfedge_const_handle()) {
if(v != Halfedge_const_handle()) {
// found shared vertex:
assert( h->vertex() == v->vertex());
assert(h->vertex() == v->vertex());
// geomtric check if the opposite segments intersect the triangles
Triangle t1( h->vertex()->point(),
Triangle t1(h->vertex()->point(),
h->next()->vertex()->point(),
h->next()->next()->vertex()->point());
Triangle t2( v->vertex()->point(),
Triangle t2(v->vertex()->point(),
v->next()->vertex()->point(),
v->next()->next()->vertex()->point());
Segment s1( h->next()->vertex()->point(),
Segment s1(h->next()->vertex()->point(),
h->next()->next()->vertex()->point());
Segment s2( v->next()->vertex()->point(),
Segment s2(v->next()->vertex()->point(),
v->next()->next()->vertex()->point());
if ( CGAL::do_intersect( t1, s2)) {
if(CGAL::do_intersect(t1, s2)) {
//cerr << "Triangles intersect (t1,s2):\n T1: " << t1
// << "\n T2 :" << t2 << endl;
triangles.push_back(t1);
triangles.push_back(t2);
} else if ( CGAL::do_intersect( t2, s1)) {
} else if(CGAL::do_intersect(t2, s1)) {
//cerr << "Triangles intersect (t2,s1):\n T1: " << t1
// << "\n T2 :" << t2 << endl;
triangles.push_back(t1);
@ -76,13 +76,13 @@ struct Intersect_facets
return;
}
// check for geometric intersection
Triangle t1( h->vertex()->point(),
Triangle t1(h->vertex()->point(),
h->next()->vertex()->point(),
h->next()->next()->vertex()->point());
Triangle t2( g->vertex()->point(),
Triangle t2(g->vertex()->point(),
g->next()->vertex()->point(),
g->next()->next()->vertex()->point());
if ( CGAL::do_intersect( t1, t2)) {
if(CGAL::do_intersect(t1, t2)) {
//cerr << "Triangles intersect:\n T1: " << t1 << "\n T2 :"
// << t2 << endl;
triangles.push_back(t1);
@ -92,26 +92,26 @@ struct Intersect_facets
};
bool Is_self_intersecting( Surface const& s )
bool Is_self_intersecting(Surface const& s)
{
std::vector<Box> boxes;
boxes.reserve( s.size_of_facets());
for ( Facet_const_iterator i = s.facets_begin(); i != s.facets_end(); ++i)
boxes.reserve(s.size_of_facets());
for(Facet_const_iterator i = s.facets_begin(); i != s.facets_end(); ++i)
{
boxes.push_back(
Box( i->halfedge()->vertex()->point().bbox()
Box(i->halfedge()->vertex()->point().bbox()
+ i->halfedge()->next()->vertex()->point().bbox()
+ i->halfedge()->next()->next()->vertex()->point().bbox(),
i));
}
std::vector<const Box*> box_ptr;
box_ptr.reserve( s.size_of_facets());
for ( std::vector<Box>::iterator j = boxes.begin(); j != boxes.end(); ++j)
box_ptr.reserve(s.size_of_facets());
for(std::vector<Box>::iterator j = boxes.begin(); j != boxes.end(); ++j)
{
box_ptr.push_back( &*j);
box_ptr.push_back(&*j);
}
CGAL::box_self_intersection_d( box_ptr.begin(), box_ptr.end(),
CGAL::box_self_intersection_d(box_ptr.begin(), box_ptr.end(),
Intersect_facets(), std::ptrdiff_t(2000));
return triangles.size() > 0 ;
return triangles.size() > 0;
}