phynexis-deps/cgal
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #1101 from maxGimeno/Enhance_triangulate_facets-GF

Browse Source 
Robustify triangulate facets
This commit is contained in:
Sebastien Loriot
2016-06-08 09:50:48 +02:00
parent d3892fe40c 1741739449
commit 1ebaebc6aa
3 changed files with 185 additions and 122 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/triangulate_faces.h
+178 -121
View File
@@ -28,7 +28,6 @@
#include <CGAL/Triangulation_vertex_base_with_info_2.h>
#include <CGAL/Triangulation_face_base_with_info_2.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <CGAL/Constrained_triangulation_plus_2.h>
#include <CGAL/Triangulation_2_projection_traits_3.h>
#include <CGAL/Polygon_mesh_processing/compute_normal.h>
@@ -75,12 +74,12 @@ class Triangulate_modifier
P_traits> Fb1;
typedef CGAL::Constrained_triangulation_face_base_2<P_traits, Fb1> Fb;
typedef CGAL::Triangulation_data_structure_2<Vb,Fb> TDS;
typedef CGAL::No_intersection_tag Itag;
typedef CGAL::Exact_intersections_tag Itag;
typedef CGAL::Constrained_Delaunay_triangulation_2<P_traits,
TDS,
Itag> CDTbase;
typedef CGAL::Constrained_triangulation_plus_2<CDTbase> CDT;
Itag> CDT;
typedef typename boost::property_traits<VertexPointMap>::reference Point_ref;
VertexPointMap _vpmap;
public:
@@ -93,126 +92,176 @@ public:
return fh->info().is_external;
}
void triangulate_face(face_descriptor f, PM& pmesh)
bool triangulate_face(face_descriptor f, PM& pmesh)
{
typename Traits::Vector_3 normal =
Polygon_mesh_processing::compute_face_normal(f, pmesh);
P_traits cdt_traits(normal);
CDT cdt(cdt_traits);
// Halfedge_around_facet_circulator
typedef typename CDT::Vertex_handle Tr_Vertex_handle;
halfedge_descriptor start = halfedge(f, pmesh);
halfedge_descriptor h = start;
Tr_Vertex_handle previous, first;
do
if(normal == typename Traits::Vector_3(0,0,0))
return false;
std::size_t original_size = CGAL::halfedges_around_face(halfedge(f, pmesh), pmesh).size();
if(original_size == 4)
{
Tr_Vertex_handle vh = cdt.insert(_vpmap[target(h, pmesh)]);
if (first == Tr_Vertex_handle()) {
first = vh;
}
vh->info() = h;
if(previous != Tr_Vertex_handle() && previous != vh) {
cdt.insert_constraint(previous, vh);
}
previous = vh;
h = next(h, pmesh);
halfedge_descriptor v0, v1, v2, v3;
v0 = halfedge(f, pmesh);
Point_ref p0 = _vpmap[target(v0, pmesh)];
v1 = next(v0, pmesh);
Point_ref p1 = _vpmap[target(v1, pmesh)];
v2 = next(v1, pmesh);
Point_ref p2 = _vpmap[target(v2, pmesh)];
v3 = next(v2, pmesh);
Point_ref p3 = _vpmap[target(v3, pmesh)];
} while( h != start );
cdt.insert_constraint(previous, first);
/* Chooses the diagonal that will split the quad in two triangles that maximize
* the scalar product of of the un-normalized normals of the two triangles.
* The lengths of the un-normalized normals (computed using cross-products of two vectors)
* are proportional to the area of the triangles.
* Maximize the scalar product of the two normals will avoid skinny triangles,
* and will also taken into account the cosine of the angle between the two normals.
* In particular, if the two triangles are oriented in different directions,
* the scalar product will be negative.
*/
double p1p3= CGAL::cross_product(p2-p1,p3-p2) * CGAL::cross_product(p0-p3,p1-p0);
double p0p2= CGAL::cross_product(p1-p0,p1-p2) * CGAL::cross_product(p3-p2,p3-p0);
// sets mark is_external
for(typename CDT::All_faces_iterator fit = cdt.all_faces_begin(),
end = cdt.all_faces_end();
fit != end; ++fit)
{
fit->info().is_external = false;
}
std::queue<typename CDT::Face_handle> face_queue;
face_queue.push(cdt.infinite_vertex()->face());
while(! face_queue.empty() )
{
typename CDT::Face_handle fh = face_queue.front();
face_queue.pop();
if(fh->info().is_external)
continue;
fh->info().is_external = true;
for(int i = 0; i <3; ++i)
if(p0p2>p1p3)
{
if(!cdt.is_constrained(typename CDT::Edge(fh, i)))
CGAL::Euler::split_face(v0, v2, pmesh);
}
else
{
CGAL::Euler::split_face(v1, v3, pmesh);
}
}
else
{
P_traits cdt_traits(normal);
CDT cdt(cdt_traits);
// Halfedge_around_facet_circulator
typedef typename CDT::Vertex_handle Tr_Vertex_handle;
halfedge_descriptor start = halfedge(f, pmesh);
halfedge_descriptor h = start;
Tr_Vertex_handle previous, first;
do
{
Tr_Vertex_handle vh = cdt.insert(_vpmap[target(h, pmesh)]);
if (first == Tr_Vertex_handle()) {
first = vh;
}
vh->info() = h;
if(previous != Tr_Vertex_handle() && previous != vh) {
cdt.insert_constraint(previous, vh);
}
previous = vh;
h = next(h, pmesh);
} while( h != start );
cdt.insert_constraint(previous, first);
// sets mark is_external
for(typename CDT::All_faces_iterator fit = cdt.all_faces_begin(),
end = cdt.all_faces_end();
fit != end; ++fit)
{
fit->info().is_external = false;
}
std::queue<typename CDT::Face_handle> face_queue;
face_queue.push(cdt.infinite_vertex()->face());
while(! face_queue.empty() )
{
typename CDT::Face_handle fh = face_queue.front();
face_queue.pop();
if(fh->info().is_external)
continue;
fh->info().is_external = true;
for(int i = 0; i <3; ++i)
{
face_queue.push(fh->neighbor(i));
if(!cdt.is_constrained(typename CDT::Edge(fh, i)))
{
face_queue.push(fh->neighbor(i));
}
}
}
if(cdt.dimension() != 2 ||
cdt.number_of_vertices() != original_size)
return false;
// then modify the polyhedron
// make_hole. (see comment in function body)
this->make_hole(halfedge(f, pmesh), pmesh);
for(typename CDT::Finite_edges_iterator eit = cdt.finite_edges_begin(),
end = cdt.finite_edges_end();
eit != end; ++eit)
{
typename CDT::Face_handle fh = eit->first;
const int index = eit->second;
typename CDT::Face_handle opposite_fh = fh->neighbor(eit->second);
const int opposite_index = opposite_fh->index(fh);
const Tr_Vertex_handle va = fh->vertex(cdt. cw(index));
const Tr_Vertex_handle vb = fh->vertex(cdt.ccw(index));
if( ! (is_external(fh) && is_external(opposite_fh))//not both fh are external
&& ! cdt.is_constrained(*eit) ) //and edge is not constrained
{
// strictly internal edge
halfedge_descriptor hnew = halfedge(add_edge(pmesh), pmesh),
hnewopp = opposite(hnew, pmesh);
fh->info().e[index] = hnew;
opposite_fh->info().e[opposite_index] = hnewopp;
set_target(hnew, target(va->info(), pmesh), pmesh);
set_target(hnewopp, target(vb->info(), pmesh), pmesh);
}
if( cdt.is_constrained(*eit) ) //edge is constrained
{
if(!is_external(fh)) {
fh->info().e[index] = va->info();
}
if(!is_external(opposite_fh)) {
opposite_fh->info().e[opposite_index] = vb->info();
}
}
}
for(typename CDT::Finite_faces_iterator fit = cdt.finite_faces_begin(),
end = cdt.finite_faces_end();
fit != end; ++fit)
{
if(!is_external(fit))
{
halfedge_descriptor h0 = fit->info().e[0];
halfedge_descriptor h1 = fit->info().e[1];
halfedge_descriptor h2 = fit->info().e[2];
CGAL_assertion(h0 != halfedge_descriptor());
CGAL_assertion(h1 != halfedge_descriptor());
CGAL_assertion(h2 != halfedge_descriptor());
set_next(h0, h1, pmesh);
set_next(h1, h2, pmesh);
set_next(h2, h0, pmesh);
Euler::fill_hole(h0, pmesh);
}
}
}
// then modify the polyhedron
// make_hole. (see comment in function body)
this->make_hole(halfedge(f, pmesh), pmesh);
for(typename CDT::Finite_edges_iterator eit = cdt.finite_edges_begin(),
end = cdt.finite_edges_end();
eit != end; ++eit)
{
typename CDT::Face_handle fh = eit->first;
const int index = eit->second;
typename CDT::Face_handle opposite_fh = fh->neighbor(eit->second);
const int opposite_index = opposite_fh->index(fh);
const Tr_Vertex_handle va = fh->vertex(cdt. cw(index));
const Tr_Vertex_handle vb = fh->vertex(cdt.ccw(index));
if( ! (is_external(fh) && is_external(opposite_fh))//not both fh are external
&& ! cdt.is_constrained(*eit) ) //and edge is not constrained
{
// strictly internal edge
halfedge_descriptor hnew = halfedge(add_edge(pmesh), pmesh),
hnewopp = opposite(hnew, pmesh);
fh->info().e[index] = hnew;
opposite_fh->info().e[opposite_index] = hnewopp;
set_target(hnew, target(va->info(), pmesh), pmesh);
set_target(hnewopp, target(vb->info(), pmesh), pmesh);
}
if( cdt.is_constrained(*eit) ) //edge is constrained
{
if(!is_external(fh)) {
fh->info().e[index] = va->info();
}
if(!is_external(opposite_fh)) {
opposite_fh->info().e[opposite_index] = vb->info();
}
}
}
for(typename CDT::Finite_faces_iterator fit = cdt.finite_faces_begin(),
end = cdt.finite_faces_end();
fit != end; ++fit)
{
if(!is_external(fit))
{
halfedge_descriptor h0 = fit->info().e[0];
halfedge_descriptor h1 = fit->info().e[1];
halfedge_descriptor h2 = fit->info().e[2];
CGAL_assertion(h0 != halfedge_descriptor());
CGAL_assertion(h1 != halfedge_descriptor());
CGAL_assertion(h2 != halfedge_descriptor());
set_next(h0, h1, pmesh);
set_next(h1, h2, pmesh);
set_next(h2, h0, pmesh);
Euler::fill_hole(h0, pmesh);
}
}
return true;
}
template<typename FaceRange>
void operator()(FaceRange face_range, PM& pmesh)
bool operator()(FaceRange face_range, PM& pmesh)
{
bool result = true;
// One need to store facet handles into a vector, because the list of
// facets of the polyhedron will be modified during the loop, and
// that invalidates the range [facets_begin(), facets_end()[.
@@ -228,8 +277,10 @@ public:
// Iterates on the vector of face descriptors
BOOST_FOREACH(face_descriptor f, facets)
{
this->triangulate_face(f, pmesh);
if(!this->triangulate_face(f, pmesh))
result = false;
}
return result;
}
void make_hole(halfedge_descriptor h, PM& pmesh)
@@ -265,14 +316,16 @@ public:
* @param pmesh the polygon mesh to which the face to be triangulated belongs
* @param np optional sequence of \ref namedparameters among the ones listed below
*
*
* \cgalNamedParamsBegin
* \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `pmesh` \cgalParamEnd
* \cgalParamBegin{geom_traits} a geometric traits class instance \cgalParamEnd
* \cgalNamedParamsEnd
*
* @return `true` if the face has been triangulated.
*/
template<typename PolygonMesh, typename NamedParameters>
void triangulate_face(typename boost::graph_traits<PolygonMesh>::face_descriptor f,
bool triangulate_face(typename boost::graph_traits<PolygonMesh>::face_descriptor f,
PolygonMesh& pmesh,
const NamedParameters& np)
{
@@ -288,14 +341,14 @@ void triangulate_face(typename boost::graph_traits<PolygonMesh>::face_descriptor
typedef typename GetGeomTraits<PolygonMesh, NamedParameters>::type Kernel;
internal::Triangulate_modifier<PolygonMesh, VPMap, Kernel> modifier(vpmap);
modifier.triangulate_face(f, pmesh);
return modifier.triangulate_face(f, pmesh);
}
template<typename PolygonMesh>
void triangulate_face(typename boost::graph_traits<PolygonMesh>::face_descriptor f,
bool triangulate_face(typename boost::graph_traits<PolygonMesh>::face_descriptor f,
PolygonMesh& pmesh)
{
triangulate_face(f, pmesh, CGAL::Polygon_mesh_processing::parameters::all_default());
return triangulate_face(f, pmesh, CGAL::Polygon_mesh_processing::parameters::all_default());
}
/**
@@ -318,9 +371,11 @@ void triangulate_face(typename boost::graph_traits<PolygonMesh>::face_descriptor
* \cgalParamBegin{geom_traits} a geometric traits class instance \cgalParamEnd
* \cgalNamedParamsEnd
*
* @return `true` if all the faces have been triangulated.
* @see triangulate_face()
*/
template <typename FaceRange, typename PolygonMesh, typename NamedParameters>
void triangulate_faces(FaceRange face_range,
bool triangulate_faces(FaceRange face_range,
PolygonMesh& pmesh,
const NamedParameters& np)
{
@@ -336,13 +391,13 @@ void triangulate_faces(FaceRange face_range,
typedef typename GetGeomTraits<PolygonMesh, NamedParameters>::type Kernel;
internal::Triangulate_modifier<PolygonMesh, VPMap, Kernel> modifier(vpmap);
modifier(face_range, pmesh);
return modifier(face_range, pmesh);
}
template <typename FaceRange, typename PolygonMesh>
void triangulate_faces(FaceRange face_range, PolygonMesh& pmesh)
bool triangulate_faces(FaceRange face_range, PolygonMesh& pmesh)
{
triangulate_faces(face_range, pmesh, CGAL::Polygon_mesh_processing::parameters::all_default());
return triangulate_faces(face_range, pmesh, CGAL::Polygon_mesh_processing::parameters::all_default());
}
/**
@@ -360,18 +415,20 @@ void triangulate_faces(FaceRange face_range, PolygonMesh& pmesh)
* \cgalParamBegin{geom_traits} a geometric traits class instance \cgalParamEnd
* \cgalNamedParamsEnd
*
* @return `true` if all the faces have been triangulated.
* @see triangulate_face()
*/
template <typename PolygonMesh, typename NamedParameters>
void triangulate_faces(PolygonMesh& pmesh,
bool triangulate_faces(PolygonMesh& pmesh,
const NamedParameters& np)
{
triangulate_faces(faces(pmesh), pmesh, np);
return triangulate_faces(faces(pmesh), pmesh, np);
}
template <typename PolygonMesh>
void triangulate_faces(PolygonMesh& pmesh)
bool triangulate_faces(PolygonMesh& pmesh)
{
triangulate_faces(faces(pmesh), pmesh, CGAL::Polygon_mesh_processing::parameters::all_default());
return triangulate_faces(faces(pmesh), pmesh, CGAL::Polygon_mesh_processing::parameters::all_default());
}
} // end namespace Polygon_mesh_processing