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

1279 lines
45 KiB
C++
Raw Blame History

// Copyright (c) 2006-2007 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Laurent RINEAU
#ifndef CGAL_POLYHEDRAL_SURFACE_3_H
#define CGAL_POLYHEDRAL_SURFACE_3_H
#include <CGAL/config.h>
#ifdef CGAL_DONT_SUBMIT
#ifdef CGAL_POLYHEDRAL_SURFACE_VERBOSE_CONSTRUCTION
# define CGAL_SURFACE_MESHER_DEBUG_POLYHEDRAL_SURFACE_CONSTRUCTION
#endif
#ifndef CGAL_SURFACE_MESHER_POLYHEDRAL_SURFACE_USE_INTERSECTION_DATA_STRUCTURE
# define CGAL_SURFACE_MESHER_POLYHEDRAL_SURFACE_USE_OCTREE 1
#endif
//#include <CGAL/Polyhedron_3.h>
#include <CGAL/enriched_polyhedron.h>
#include <CGAL/HalfedgeDS_vector.h>
#include <CGAL/IO/Polyhedron_iostream.h>
#include <CGAL/Inverse_index.h>
#include <CGAL/circulator.h>
#include <CGAL/iterator.h>
#include <memory>
#include <boost/tuple/tuple.hpp>
#include <CGAL/make_surface_mesh.h>
#ifdef CGAL_SURFACE_MESHER_POLYHEDRAL_SURFACE_USE_OCTREE
# include <CGAL/Data_structure_using_octree_3.h>
#endif
#ifdef CGAL_SURFACE_MESHER_POLYHEDRAL_SURFACE_USE_INTERSECTION_DATA_STRUCTURE
# include <CGAL/Surface_mesher/Intersection_data_structure_3.h>
#endif
#ifdef CGAL_SURFACE_MESHER_POLYHEDRAL_SURFACE_USE_PINPOLYHEDRON
# include <pinpolyhedron.h>
#endif
#include <CGAL/Surface_mesher/Polyhedral_oracle.h>
#include <CGAL/Surface_mesher/Has_edges.h>
#include <iostream>
#include <vector>
#include <set>
#include <queue>
#include <algorithm> // random_shuffle, distance, copy, set_intersection
#ifdef CGAL_SURFACE_MESHER_DEBUG_POLYHEDRAL_SURFACE_CONSTRUCTION
#include <boost/format.hpp>
#include <CGAL/Timer.h>
#endif // CGAL_SURFACE_MESHER_DEBUG_POLYHEDRAL_SURFACE_CONSTRUCTION
namespace CGAL {
namespace Surface_mesher {
template <
class Surface,
class Point_creator,
class Visitor,
class Tag,
bool
>
class Polyhedral_oracle;
} // end namespace Surface_mesher
template <class Polyhedron>
struct Polyhedral_surface_3_type_helper {
typedef typename Polyhedron::Vertex_handle Vertex_handle;
typedef typename Polyhedron::Facet_handle Facet_handle;
typedef typename Polyhedron::Halfedge_handle Halfedge_handle;
typedef Halfedge_handle Edge_handle;
typedef typename Polyhedron::Traits::Point_3 Point_3;
struct Vertex_node {
typedef std::set<int> Incident_edges; // set of edge id
Incident_edges incident_edges;
};
struct Edge_node {
typedef std::vector<Edge_handle> Sub_edges;
typedef std::set<int> Incident_facets; // set of facet id
Sub_edges sub_edges;
Incident_facets incident_facets;
};
struct Facet_node {
typedef std::vector<Facet_handle> Sub_facets;
Sub_facets sub_facets;
};
struct Incidence_graph {
typedef std::vector<Vertex_node> Vertices_nodes;
typedef std::vector<Edge_node> Edges_nodes;
typedef std::vector<Facet_node> Facet_nodes;
typedef typename Vertices_nodes::iterator Vertex_iterator;
typedef typename Edges_nodes::iterator Edge_iterator;
typedef typename Facet_nodes::iterator Facet_iterator;
Vertices_nodes vertices;
Edges_nodes edges;
Facet_nodes facets;
void clear()
{
vertices.clear();
edges.clear();
facets.clear();
}
int vertex_id(Vertex_iterator vit)
{
return std::distance(vertices.begin(), vit);
}
int edge_id(Edge_iterator eit)
{
return std::distance(edges.begin(), eit);
}
int facet_id(Facet_iterator fit)
{
return std::distance(facets.begin(), fit);
}
void gl_draw_edge(const int edge_index)
{
::glBegin(GL_LINES);
typedef typename Edge_node::Sub_edges Sub_edges;
for(typename Sub_edges::iterator
eit = edges[edge_index].sub_edges.begin(),
end = edges[edge_index].sub_edges.end();
eit != end;
++eit)
{
const Point_3& p1 = (*eit)->opposite()->vertex()->point();
const Point_3& p2 = (*eit)->vertex()->point();
::glVertex3f(p1[0],p1[1],p1[2]);
::glVertex3f(p2[0],p2[1],p2[2]);
}
::glEnd();
} // end gl_draw_edge(...)
void gl_draw_facet(const int facet_index,
bool smooth_shading,
bool use_normals,
bool inverse_normals = false)
{
// draw triangles
::glBegin(GL_TRIANGLES);
typedef typename Facet_node::Sub_facets Sub_facets;
for(typename Sub_facets::iterator
fit = facets[facet_index].sub_facets.begin(),
end_fit = facets[facet_index].sub_facets.end();
fit != end_fit;
++fit)
{
// one normal per face
if(use_normals && !smooth_shading)
{
typedef typename Polyhedron::Facet Facet;
const typename Facet::Normal_3& normal = (*fit)->normal();
if(inverse_normals)
::glNormal3f(-normal[0],-normal[1],-normal[2]);
else
::glNormal3f(normal[0],normal[1],normal[2]);
}
typename Polyhedron::Halfedge_around_facet_circulator
he = (*fit)->facet_begin(),
end_he = he;
do
{
// one normal per vertex
if(use_normals && smooth_shading)
{
const typename Polyhedron::Facet::Normal_3& normal = he->vertex()->normal();
if(inverse_normals)
::glNormal3f(-normal[0],-normal[1],-normal[2]);
else
::glNormal3f(normal[0],normal[1],normal[2]);
}
// polygon assembly is performed per vertex
const Point_3& point = he->vertex()->point();
::glVertex3d(point[0],point[1],point[2]);
}
while(++he != end_he);
}
::glEnd(); // end polygon assembly
} // end gl_draw_facet(...)
}; // end Incidence_graph
}; // end Polyhedral_surface_3_type_helper
template <class GT,
class Has_edges_tag_ = Surface_mesher::Has_no_edges,
class Polyhedron_3 = Enriched_polyhedron<GT,
Enriched_items> >
class Polyhedral_surface_3 :
public Has_edges_tag_,
public Polyhedron_3
{
public:
typedef GT Geom_traits;
typedef Has_edges_tag_ Has_edges_tag;
typedef Polyhedron_3 Polyhedron;
typedef Polyhedral_surface_3<Geom_traits,
Has_edges_tag,
Polyhedron> Self;
typedef typename Polyhedron::Vertex_handle Vertex_handle;
typedef typename Polyhedron::Facet_handle Facet_handle;
typedef typename Polyhedron::Halfedge_handle Halfedge_handle;
typedef Halfedge_handle Edge_handle;
typedef typename Polyhedron::Halfedge_around_vertex_circulator
Halfedge_around_vertex_circulator;
typedef typename Polyhedron::Halfedge_around_facet_circulator
Halfedge_around_facet_circulator;
typedef typename Polyhedron::Halfedge_iterator Halfedge_iterator;
typedef typename Polyhedron::Facet_iterator Facet_iterator;
struct Compare_vertex_iterators {
bool operator()(const typename Polyhedron_3::Vertex_iterator& va,
const typename Polyhedron_3::Vertex_iterator& vb) const
{
return (&*va)<(&*vb);
}
};
using Polyhedron::halfedges_begin;
using Polyhedron::halfedges_end;
using Polyhedron::edges_begin;
using Polyhedron::edges_end;
using Polyhedron::facets_begin;
using Polyhedron::facets_end;
typedef Polyhedral_surface_3_type_helper<Polyhedron> Type_helper;
typedef typename Type_helper::Vertex_node Graph_vertex_node;
typedef typename Type_helper::Edge_node Graph_edge_node;
typedef typename Type_helper::Facet_node Graph_facet_node;
typedef typename Type_helper::Incidence_graph Incidence_graph;
class Normalized_geom_traits : public Geom_traits
{
public:
typedef typename
Kernel_traits<typename Geom_traits::Point_3>::Kernel::Point_3 Point_3;
};
typedef typename Geom_traits::FT FT;
typedef typename Geom_traits::Point_3 Point_3;
typedef typename Geom_traits::Segment_3 Segment_3;
typedef typename Geom_traits::Triangle_3 Triangle_3;
typedef typename Geom_traits::Vector_3 Vector_3;
typedef Bbox_3 Bbox;
#ifdef CGAL_SURFACE_MESHER_POLYHEDRAL_SURFACE_USE_OCTREE
typedef Data_structure_using_octree_3<Normalized_geom_traits> Subfacets_tree;
typedef Data_structure_using_octree_3<Normalized_geom_traits> Subsegments_tree;
#endif
#ifdef CGAL_SURFACE_MESHER_POLYHEDRAL_SURFACE_USE_INTERSECTION_DATA_STRUCTURE
typedef Intersection_data_structure_3<Normalized_geom_traits,
Triangle_3> Subfacets_tree;
typedef Intersection_data_structure_3<Normalized_geom_traits,
Segment_3> Subsegments_tree;
#endif
typedef std::shared_ptr<Subfacets_tree> Subfacets_tree_ptr;
typedef std::shared_ptr<Subsegments_tree> Subsegments_tree_ptr;
// typedef typename Subsegments_tree::Point_with_index
// Intersection_point;
template <
class Surface,
class Point_creator,
class Visitor,
class Tag
>
friend class Surface_mesher::Polyhedral_oracle;
typedef Surface_mesher::Polyhedral_oracle<Self> Surface_mesher_traits_3;
Polyhedral_surface_3(const double sharp_edges_angle_lower_bound = 60.,
const double sharp_edges_angle_upper_bound = 180.,
const double sharp_vertices_angle_lower_bound = 30.,
const double sharp_vertices_angle_upper_bound = 180.) :
subfacets_tree_ptr(new Subfacets_tree()),
#ifdef CGAL_SURFACE_MESHER_POLYHEDRAL_SURFACE_USE_OCTREE
subsegments_tree_ptr(new Subsegments_tree(false, true, false)),
#endif
#ifdef CGAL_SURFACE_MESHER_POLYHEDRAL_SURFACE_USE_INTERSECTION_DATA_STRUCTURE
subsegments_tree_ptr(new Subsegments_tree()),
#endif
input_vertices_ptr(new Input_vertices()),
corner_vertices_ptr(new Corner_vertices()),
edges_vertices_ptr(new Edges_vertices()),
sharp_edges_angle_lower_bound(sharp_edges_angle_lower_bound),
sharp_edges_angle_upper_bound(sharp_edges_angle_upper_bound),
sharp_vertices_angle_lower_bound(sharp_vertices_angle_lower_bound),
sharp_vertices_angle_upper_bound(sharp_vertices_angle_upper_bound)
{
}
Polyhedral_surface_3(std::istream& input_file,
const double sharp_edges_angle_lower_bound = 60.,
const double sharp_edges_angle_upper_bound = 180.,
bool auto_construct_octree = true) :
subfacets_tree_ptr(new Subfacets_tree()),
#ifdef CGAL_SURFACE_MESHER_POLYHEDRAL_SURFACE_USE_OCTREE
subsegments_tree_ptr(new Subsegments_tree(false, true, false)),
#endif
#ifdef CGAL_SURFACE_MESHER_POLYHEDRAL_SURFACE_USE_INTERSECTION_DATA_STRUCTURE
subsegments_tree_ptr(new Subsegments_tree()),
#endif
input_vertices_ptr(new Input_vertices()),
corner_vertices_ptr(new Corner_vertices()),
edges_vertices_ptr(new Edges_vertices()),
sharp_edges_angle_lower_bound(sharp_edges_angle_lower_bound),
sharp_edges_angle_upper_bound(sharp_edges_angle_upper_bound)
{
#ifdef CGAL_SURFACE_MESHER_DEBUG_POLYHEDRAL_SURFACE_CONSTRUCTION
CGAL::Timer timer;
std::cerr << "Creating polyhedron... ";
timer.start();
#endif // CGAL_SURFACE_MESHER_DEBUG_POLYHEDRAL_SURFACE_CONSTRUCTION
CGAL::scan_OFF(input_file, *this, true);
if(!input_file)
return;
this->compute_bounding_box();
this->compute_normals();
bounding_box = Bbox(this->xmin(),
this->ymin(),
this->zmin(),
this->xmax(),
this->ymax(),
this->zmax());
bounding_box_sq_radius = bounding_box.xmax()-bounding_box.xmin();
bounding_box_sq_radius =
CGAL_NTS max BOOST_PREVENT_MACRO_SUBSTITUTION
(bounding_box_sq_radius,
FT(bounding_box.ymax()-bounding_box.ymin()));
bounding_box_sq_radius =
CGAL_NTS max BOOST_PREVENT_MACRO_SUBSTITUTION
(bounding_box_sq_radius,
FT(bounding_box.zmax()-bounding_box.zmin()));
bounding_box_sq_radius /= 2;
bounding_box_sq_radius *= bounding_box_sq_radius;
bounding_box_sq_radius *= 3;
#ifdef CGAL_SURFACE_MESHER_DEBUG_POLYHEDRAL_SURFACE_CONSTRUCTION
timer.stop();
std::cerr <<
::boost::format("done (%1%s)\n"
" number of vertices: %2%\n"
" number of edges: %3%\n"
" number of facets: %4%\n")
% timer.time()
% this->size_of_vertices()
% ( this->size_of_halfedges() / 2 )
% this->size_of_facets();
#endif // CGAL_SURFACE_MESHER_DEBUG_POLYHEDRAL_SURFACE_CONSTRUCTION
if(auto_construct_octree)
{
compute_sharp_edges_incidence_graph();
construct_octree();
}
} // end of Polyhedral_surface_3 constructor
void construct_octree()
{
#if 0
GT gt = GT();
typename GT::Construct_orthogonal_vector_3 orthogonal_vector =
gt.construct_orthogonal_vector_3_object();
typename GT::Compute_squared_length_3 squared_length =
gt.compute_squared_length_3_object();
typename GT::Compute_scalar_product_3 scalar_product =
gt.compute_scalar_product_3_object();
// typedef CGAL::Polyhedron_3<GT> Polyhedron_3;
// typedef Enriched_polyhedron<GT,Enriched_items> Polyhedron_3;
class Facet_ortho_vector {
GT gt;
public:
Facet_ortho_vector(GT gt) : gt(gt) {}
typename GT::Vector_3
operator()(const typename Polyhedron_3::Facet& f) {
typename Polyhedron_3::Halfedge_around_facet_const_circulator
edges_circ = f.facet_begin();
const Point_3& p1 = edges_circ++->vertex()->point();
const Point_3& p2 = edges_circ++->vertex()->point();
const Point_3& p3 = edges_circ++->vertex()->point();
return gt.construct_orthogonal_vector_3_object()(p1, p2, p3);
}
};
Facet_ortho_vector facet_ortho_vector(gt);
const double cosine_bound =
std::cos(CGAL_PI * sharp_edges_angle_lower_bound / 360. );
const FT cosine_squared_bound = cosine_bound * cosine_bound;
#endif // if 0
typedef typename Polyhedron_3::Vertex_const_iterator Vertex_const_it;
#ifdef CGAL_SURFACE_MESHER_POLYHEDRAL_SURFACE_USE_PINPOLYHEDRON
double (*vertices_array)[3] =
(double (*)[3])new double [3*this->size_of_vertices()];
// static_cast<double(*)[3]>(new double [3*this->size_of_vertices()]);
int (*facets_array)[3] =
(int (*)[3])new int [3*this->size_of_facets()];
int i = 0;
typedef Inverse_index<Vertex_const_it> Index;
Index index( this->vertices_begin(), this->vertices_end());
#endif // CGAL_SURFACE_MESHER_POLYHEDRAL_SURFACE_USE_PINPOLYHEDRON
for(typename Polyhedron_3::Vertex_const_iterator
vit = this->vertices_begin(),
end = this->vertices_end();
vit != end;
++vit)
{
#ifdef CGAL_SURFACE_MESHER_POLYHEDRAL_SURFACE_USE_OCTREE
subfacets_tree_ptr->add_constrained_vertex(vit->point());
#endif
// input_vertices_ptr->insert(vit);
#ifdef CGAL_SURFACE_MESHER_POLYHEDRAL_SURFACE_USE_PINPOLYHEDRON
vertices_array[i][0] = vit->point().x();
vertices_array[i][1] = vit->point().y();
vertices_array[i][2] = vit->point().z();
++i;
#endif // CGAL_SURFACE_MESHER_POLYHEDRAL_SURFACE_USE_PINPOLYHEDRON
}
typename Polyhedron_3::size_type facet_index = 0;
for(typename Polyhedron_3::Facet_const_iterator fit =
this->facets_begin();
fit != this->facets_end();
++fit, ++facet_index)
{
CGAL_assertion(fit->is_triangle());
typename Polyhedron_3::Halfedge_around_facet_const_circulator
edges_circ = fit->facet_begin();
const Point_3& p1 = edges_circ++->vertex()->point();
const Point_3& p2 = edges_circ++->vertex()->point();
const Point_3& p3 = edges_circ++->vertex()->point();
#ifdef CGAL_SURFACE_MESHER_POLYHEDRAL_SURFACE_USE_PINPOLYHEDRON
edges_circ = fit->facet_begin();
facets_array[facet_index][0] = index[edges_circ++->vertex()];
facets_array[facet_index][1] = index[edges_circ++->vertex()];
facets_array[facet_index][2] = index[edges_circ++->vertex()];
#endif // CGAL_SURFACE_MESHER_POLYHEDRAL_SURFACE_USE_PINPOLYHEDRON
#ifdef CGAL_SURFACE_MESHER_POLYHEDRAL_SURFACE_USE_OCTREE
CGAL_assertion(fit->tag() >= 0);
subfacets_tree_ptr->add_constrained_facet(p1, p2, p3, fit->tag());
#endif
#ifdef CGAL_SURFACE_MESHER_POLYHEDRAL_SURFACE_USE_INTERSECTION_DATA_STRUCTURE
subfacets_tree_ptr->add_element(Triangle_3(p1, p2, p3));
#endif
#ifdef CGAL_POLYHEDRAL_SURFACE_VERBOSE_CONSTRUCTION
std::cerr << ::boost::format("new facet: #%4% (%1%, %2%, %3%)\n")
% p1 % p2 % p3 % facet_index;
#endif // CGAL_POLYHEDRAL_SURFACE_VERBOSE_CONSTRUCTION
}
#ifdef CGAL_SURFACE_MESHER_POLYHEDRAL_SURFACE_USE_PINPOLYHEDRON
# ifdef CGAL_SURFACE_MESHER_DEBUG_POLYHEDRAL_SURFACE_CONSTRUCTION
std::cerr << "Creating the PointInPolyhedron data structure... ";
timer.reset();
timer.start();
# endif // CGAL_SURFACE_MESHER_DEBUG_POLYHEDRAL_SURFACE_CONSTRUCTION
pinpolyhedron_ptr =
PointInPolyhedron_ptr(new PointInPolyhedron(vertices_array,
this->size_of_vertices(),
facets_array,
this->size_of_facets()));
delete [] vertices_array;
delete [] facets_array;
# ifdef CGAL_SURFACE_MESHER_DEBUG_POLYHEDRAL_SURFACE_CONSTRUCTION
timer.stop();
std::cerr << ::boost::format("done (%1%s)\n") % timer.time();
# endif // CGAL_SURFACE_MESHER_DEBUG_POLYHEDRAL_SURFACE_CONSTRUCTION
#endif // CGAL_SURFACE_MESHER_POLYHEDRAL_SURFACE_USE_PINPOLYHEDRON
#ifdef CGAL_SURFACE_MESHER_DEBUG_POLYHEDRAL_SURFACE_CONSTRUCTION
std::cerr << "Creating subfacets_tree... ";
CGAL::Timer timer;
timer.start();
#endif // CGAL_SURFACE_MESHER_DEBUG_POLYHEDRAL_SURFACE_CONSTRUCTION
subfacets_tree_ptr->create_data_structure();
#ifdef CGAL_SURFACE_MESHER_DEBUG_POLYHEDRAL_SURFACE_CONSTRUCTION
timer.stop();
# ifdef CGAL_SURFACE_MESHER_POLYHEDRAL_SURFACE_USE_OCTREE
std::cerr <<
::boost::format("done (%1%s)\n"
" number of facets: %3%\n"
" number of constraints: %4% (in subfacets_tree)\n")
% timer.time()
% subfacets_tree_ptr->number_of_vertices()
% subfacets_tree_ptr->number_of_facets()
% subfacets_tree_ptr->number_of_constraints();
# endif // CGAL_SURFACE_MESHER_POLYHEDRAL_SURFACE_USE_OCTREE
# ifdef CGAL_SURFACE_MESHER_POLYHEDRAL_SURFACE_USE_INTERSECTION_DATA_STRUCTURE
std::cerr <<
::boost::format("done (%1%s)\n"
" number of facets: %2%\n")
% timer.time()
% subfacets_tree_ptr->number_of_elements();
# endif // CGAL_SURFACE_MESHER_POLYHEDRAL_SURFACE_USE_INTERSECTION_DATA_STRUCTURE
#endif // CGAL_SURFACE_MESHER_DEBUG_POLYHEDRAL_SURFACE_CONSTRUCTION
if( this->has_edges() ) {
for(typename Polyhedron_3::Edge_const_iterator eit =
this->edges_begin();
eit != this->edges_end();
++eit)
{
if(eit->tag() >= 0)
{
const Point_3 pa = eit->vertex()->point();
const Point_3 pb = eit->opposite()->vertex()->point();
#ifdef CGAL_SURFACE_MESHER_POLYHEDRAL_SURFACE_USE_OCTREE
subsegments_tree_ptr->add_constrained_edge(pa, pb, eit->tag());
#endif
#ifdef CGAL_SURFACE_MESHER_POLYHEDRAL_SURFACE_USE_INTERSECTION_DATA_STRUCTURE
subsegments_tree_ptr->add_element(Segment_3(pa, pb));
# ifdef CGAL_POLYHEDRAL_SURFACE_VERBOSE_CONSTRUCTION
std::cerr <<
::boost::format("new edge: (%1%, %2%)")
% pa % pb;
if(eit->is_border_edge())
std:: cerr << " (on border)\n";
else
{
struct Triangle {
typename GT::Triangle_3 operator()(typename Polyhedron_3::Facet facet) {
CGAL_assertion(facet.is_triangle());
typename Polyhedron_3::Halfedge_around_facet_const_circulator
edges_circ = facet.facet_begin();
const Point_3& p1 = edges_circ++->vertex()->point();
const Point_3& p2 = edges_circ++->vertex()->point();
const Point_3& p3 = edges_circ++->vertex()->point();
return typename GT::Triangle_3(p1, p2, p3);
}
};
std:: cerr << " (not on border) ";
std::cerr <<
::boost::format("facets=(%1%, %2%)")
% Triangle()(*(eit->facet()))
% Triangle()(*(eit->opposite()->facet()));
}
# endif // CGAL_POLYHEDRAL_SURFACE_VERBOSE_CONSTRUCTION
#endif // CGAL_SURFACE_MESHER_POLYHEDRAL_SURFACE_USE_OCTREE
}
}
#if 0
typedef std::map<Point_3, int> Edges_vertex_counter;
Edges_vertex_counter edges_vertex_counter;
for(typename Polyhedron_3::Edge_const_iterator eit =
this->edges_begin();
eit != this->edges_end();
++eit)
{
typename Polyhedron_3::Halfedge_const_handle opposite = eit->opposite();
bool insert_that_edge = false;
if( eit->is_border_edge() )
insert_that_edge = true;
else
{
// CGAL_assertion(eit->is_triangle());
// CGAL_assertion(opposite->is_triangle());
if(eit->facet()->facet_degree() != 3)
std::cerr << "warning: degree=" << eit->facet_degree() << "\n";
if(opposite->facet()->facet_degree() != 3)
std::cerr << "warning: degree(opposite)=" << opposite->facet_degree() << "\n";
const Vector_3 v1 = facet_ortho_vector(*eit->facet());
const Vector_3 v2 = facet_ortho_vector(*eit->opposite()->facet());
const FT product = scalar_product(v1, v2);
insert_that_edge = product < FT(0) ||
product * product <
cosine_squared_bound * squared_length(v1) * squared_length(v2);
}
if(insert_that_edge)
{
const Point_3 pa = eit->vertex()->point();
const Point_3 pb = opposite->vertex()->point();
++edges_vertex_counter[pa];
++edges_vertex_counter[pb];
#ifdef CGAL_SURFACE_MESHER_POLYHEDRAL_SURFACE_USE_OCTREE
subsegments_tree_ptr->add_constrained_edge(pa, pb);
#endif
#ifdef CGAL_SURFACE_MESHER_POLYHEDRAL_SURFACE_USE_INTERSECTION_DATA_STRUCTURE
subsegments_tree_ptr->add_element(Segment_3(pa, pb));
# ifdef CGAL_POLYHEDRAL_SURFACE_VERBOSE_CONSTRUCTION
std::cerr <<
::boost::format("new edge: (%1%, %2%)")
% pa % pb;
if(eit->is_border_edge())
std:: cerr << " (on border)\n";
else
{
struct Triangle {
typename GT::Triangle_3 operator()(typename Polyhedron_3::Facet facet) {
CGAL_assertion(facet.is_triangle());
typename Polyhedron_3::Halfedge_around_facet_const_circulator
edges_circ = facet.facet_begin();
const Point_3& p1 = edges_circ++->vertex()->point();
const Point_3& p2 = edges_circ++->vertex()->point();
const Point_3& p3 = edges_circ++->vertex()->point();
return typename GT::Triangle_3(p1, p2, p3);
}
};
std:: cerr << " (not on border) ";
std::cerr <<
::boost::format("facets=(%3%, %4%) normals=(%1%, %2%)\n")
% facet_ortho_vector(*eit->facet())
% facet_ortho_vector(*eit->opposite()->facet())
% Triangle()(*(eit->facet()))
% Triangle()(*(eit->opposite()->facet()));
}
# endif // CGAL_POLYHEDRAL_SURFACE_VERBOSE_CONSTRUCTION
#endif // CGAL_SURFACE_MESHER_POLYHEDRAL_SURFACE_USE_OCTREE
}
}
for(typename Edges_vertex_counter::const_iterator it =
edges_vertex_counter.begin();
it != edges_vertex_counter.end();
++it)
{
input_vertices_ptr->erase(it->first);
if(it->second != 2)
{
#ifdef CGAL_POLYHEDRAL_SURFACE_VERBOSE_CONSTRUCTION
std::cerr << ::boost::format("corner point: (%1%)\n")
% it->first;
#endif // CGAL_POLYHEDRAL_SURFACE_VERBOSE_CONSTRUCTION
corner_vertices_ptr->push_back(it->first);
}
else
{
#ifdef CGAL_POLYHEDRAL_SURFACE_VERBOSE_CONSTRUCTION
std::cerr << ::boost::format("edge point: (%1%)\n")
% it->first;
#endif // CGAL_POLYHEDRAL_SURFACE_VERBOSE_CONSTRUCTION
edges_vertices_ptr->push_back(it->first);
}
}
// if(!corner_vertices_ptr->empty() && edges_vertices_ptr->empty())
// {
// #ifdef CGAL_SURFACE_MESHER_DEBUG_POLYHEDRAL_SURFACE_CONSTRUCTION
// std::cerr << "Incorrect input data. "
// << "Swap corner vertices and edges vertices...\n";
// #endif
// std::swap(corner_vertices, edges_vertices);
// }
#ifdef CGAL_SURFACE_MESHER_DEBUG_POLYHEDRAL_SURFACE_CONSTRUCTION
std::cerr << "Shuffle edges vertices... ";
#endif
std::random_shuffle(edges_vertices_ptr->begin(), edges_vertices_ptr->end());
#ifdef CGAL_SURFACE_MESHER_DEBUG_POLYHEDRAL_SURFACE_CONSTRUCTION
std::cerr << "done\n";
#endif
#ifdef CGAL_SURFACE_MESHER_DEBUG_POLYHEDRAL_SURFACE_CONSTRUCTION
std::cerr <<
::boost::format("number of corner vertices: %1%\n"
"number of edges vertices: %2%\n")
% corner_vertices_ptr->size()
% edges_vertices_ptr->size();
#endif // CGAL_SURFACE_MESHER_DEBUG_POLYHEDRAL_SURFACE_CONSTRUCTION
#endif // if 0
#ifdef CGAL_SURFACE_MESHER_DEBUG_POLYHEDRAL_SURFACE_CONSTRUCTION
std::cerr << "Creating subsegments_tree... ";
timer.reset();
timer.start();
#endif // CGAL_SURFACE_MESHER_DEBUG_POLYHEDRAL_SURFACE_CONSTRUCTION
subsegments_tree_ptr->create_data_structure();
#ifdef CGAL_SURFACE_MESHER_DEBUG_POLYHEDRAL_SURFACE_CONSTRUCTION
timer.stop();
# ifdef CGAL_SURFACE_MESHER_POLYHEDRAL_SURFACE_USE_OCTREE
std::cerr <<
::boost::format("done (%1%s)\n"
" number of edges: %3%\n"
" number of constraints: %4% (in subsegments_tree)\n")
% timer.time()
% subsegments_tree_ptr->number_of_vertices()
% subsegments_tree_ptr->number_of_edges()
% subsegments_tree_ptr->number_of_constraints();
# endif // CGAL_SURFACE_MESHER_POLYHEDRAL_SURFACE_USE_OCTREE
# ifdef CGAL_SURFACE_MESHER_POLYHEDRAL_SURFACE_USE_INTERSECTION_DATA_STRUCTURE
std::cerr <<
::boost::format("done (%1%s)\n"
" number of edges: %2%\n")
% timer.time()
% subsegments_tree_ptr->number_of_elements();
# endif // CGAL_SURFACE_MESHER_POLYHEDRAL_SURFACE_USE_INTERSECTION_DATA_STRUCTURE
#endif // CGAL_SURFACE_MESHER_DEBUG_POLYHEDRAL_SURFACE_CONSTRUCTION
} // end "if(this->has_edges())"
// subfacets_tree_ptr->input(input_file,
// std::back_inserter(input_vertices));
} // end construct_octree()
void set_sharp_edges_angle_bounds(double lower_bound,
double upper_bound) {
sharp_edges_angle_lower_bound = lower_bound;
sharp_edges_angle_upper_bound = upper_bound;
}
void set_sharp_vertices_angle_bounds(double lower_bound,
double upper_bound) {
sharp_vertices_angle_lower_bound = lower_bound;
sharp_vertices_angle_upper_bound = upper_bound;
}
unsigned int tag_border_edges()
{
unsigned int nb = 0;
for(Halfedge_iterator he = edges_begin();
he != edges_end();
he++)
{
const bool tag = ( he->sharp() ||
he->is_border() ||
he->opposite()->is_border() );
he->sharp() = tag;
he->opposite()->sharp() = tag;
nb += tag ? 1 : 0;
}
return nb;
}
void compute_sharp_edges_incidence_graph()
{
this->tag_sharp_edges(sharp_edges_angle_lower_bound);
this->tag_border_edges();
construct_incidence_graph();
}
const Bbox& bbox() const
{
return bounding_box;
}
const FT& bounding_sphere_squared_radius() const
{
return bounding_box_sq_radius;
}
static void new_sub_edge(Graph_edge_node& edge_node,
const int edge_index,
Halfedge_handle& he)
{
he->tag(edge_index);
he->opposite()->tag(edge_index);
edge_node.sub_edges.push_back(he);
}
static void new_sub_facet(Graph_facet_node& facet_node,
const int facet_index,
Facet_handle& fh)
{
fh->tag(facet_index);
facet_node.sub_facets.push_back(fh);
}
Vertex_handle follow_the_edge(Graph_edge_node& edge_node,
const int edge_index,
Halfedge_handle he) const
{
using boost::tie;
bool can_continue;
Halfedge_handle next_he;
tie(can_continue, next_he) = can_follow_sharp_edges(he);
while(can_continue && next_he->tag() < 0)
{
CGAL_assertion(he->vertex()->tag()<0);
he->vertex()->tag(edge_index);
new_sub_edge(edge_node, edge_index, next_he);
he = next_he;
tie(can_continue, next_he) = Self::can_follow_sharp_edges(he);
}
if(can_continue)
{
he->vertex()->tag(edge_index);
return Vertex_handle();
}
else
return he->vertex();
}
void setup_incident_vertex_to_an_edge(const int edge_index,
const Vertex_handle& v)
{
if(v->tag() < 0)
{
Graph_vertex_node vertex_node;
vertex_node.incident_edges.insert(edge_index);
v->tag(incidence_graph.vertices.size());
incidence_graph.vertices.push_back(vertex_node);
}
else
incidence_graph.vertices[v->tag()].incident_edges.insert(edge_index);
}
void construct_incidence_graph()
{
incidence_graph.clear();
input_vertices_ptr->clear();
corner_vertices_ptr->clear();
edges_vertices_ptr->clear();
this->tag_vertices(-1);
this->tag_halfedges(-1);
this->tag_facets(-1);
for(Halfedge_iterator he = halfedges_begin();
he != halfedges_end();
he++)
{
if(he->sharp() && he->tag() < 0)
{
const int edge_index = incidence_graph.edges.size();
Graph_edge_node edge_node;
new_sub_edge(edge_node, edge_index, he);
const Vertex_handle& v1 = follow_the_edge(edge_node,
edge_index,
he);
if(v1 != Vertex_handle())
setup_incident_vertex_to_an_edge(edge_index, v1);
const Vertex_handle& v2 = follow_the_edge(edge_node,
edge_index,
he->opposite());
if(v2 != Vertex_handle())
setup_incident_vertex_to_an_edge(edge_index, v2);
incidence_graph.edges.push_back(edge_node);
} // end "if halfedge he is sharp and not tagged"
} // end "for all halfedges"
for(Facet_iterator fit = facets_begin(), end = facets_end();
fit != end;
++fit)
{
if(fit->tag() < 0)
{
const int facet_index = incidence_graph.facets.size();
Graph_facet_node facet_node;
construct_graph_facet(fit,
facet_index,
facet_node,
incidence_graph);
incidence_graph.facets.push_back(facet_node);
}
}
// display to std::cerr, for debugging
std::cerr << "INCIDENCE GRAPH SUMMARY\n";
std::cerr << "Number of vertices: " << incidence_graph.vertices.size()
<< "\nNumber of edges: " << incidence_graph.edges.size()
<< "\nNumber of facets: " << incidence_graph.facets.size()
<< std::endl;
std::cerr << "VERTICES\n";
for(int i = 0, end_i = incidence_graph.vertices.size(); i < end_i; ++i)
{
std::cerr << boost::format("Vertex #%1% incident edges: ") % i;
for(typename Graph_vertex_node::Incident_edges::iterator
sub_edge_id_it = incidence_graph.vertices[i].incident_edges.begin(),
end = incidence_graph.vertices[i].incident_edges.end();
sub_edge_id_it != end; ++sub_edge_id_it)
{
std::cerr <<
boost::format("%1% ") % *sub_edge_id_it;
}
std::cerr << std::endl;
}
std::cerr << "EDGES\n";
for(int i = 0, end_i = incidence_graph.edges.size(); i < end_i; ++i)
{
std::cerr <<
boost::format("Edge #%1% number of sub-edges: %2%\n")
% i % incidence_graph.edges[i].sub_edges.size();
std::cerr <<
boost::format("Edge #%1% incident_facets: ") % i;
for(typename Graph_edge_node::Incident_facets::iterator
facet_id_it = incidence_graph.edges[i].incident_facets.begin(),
end = incidence_graph.edges[i].incident_facets.end();
facet_id_it != end; ++facet_id_it)
{
std::cerr <<
boost::format("%1% ") % *facet_id_it;
}
std::cerr << std::endl;
}
std::cerr << "FACETS\n";
for(int i = 0, end_i = incidence_graph.facets.size(); i < end_i; ++i)
{
std::cerr <<
boost::format("Facet #%1% number of sub-facets: %2%\n")
% i % incidence_graph.facets[i].sub_facets.size();
}
for(typename Polyhedron_3::Vertex_iterator
vit = this->vertices_begin(),
end = this->vertices_end();
vit != end;
++vit)
{
switch(this->type(vit)) {
case Polyhedron::SMOOTH:
vit->tag(vit->halfedge()->facet()->tag());
CGAL_assertion(vit->tag()>=0);
input_vertices_ptr->insert(vit);
break;
case Polyhedron::DART:
case Polyhedron::CORNER:
CGAL_assertion(vit->tag()>=0);
corner_vertices_ptr->push_back(vit);
case Polyhedron::CREASE_REGULAR:
case Polyhedron::CREASE_IRREGULAR:
edges_vertices_ptr->push_back(vit);
}
}
}
void construct_graph_facet(Facet_handle seed_facet,
const int facet_index,
Graph_facet_node& facet_node,
Incidence_graph& graph)
{
typedef std::queue<Facet_handle> Facets_queue;
Facets_queue facets_queue;
facets_queue.push(seed_facet);
while(!facets_queue.empty())
{
Facet_handle fh = facets_queue.front();
facets_queue.pop();
if(fh->tag() < 0)
{
fh->tag(facet_index);
facet_node.sub_facets.push_back(fh);
Halfedge_around_facet_circulator fit = fh->facet_begin();
Halfedge_around_facet_circulator end = fit;
CGAL_For_all(fit,end)
{
const int edge_index = fit->tag();
const Halfedge_handle& neighbor_facet_he = fit->opposite();
if(!neighbor_facet_he->is_border())
{
if(edge_index < 0) // if not yet handled...
{
facets_queue.push(neighbor_facet_he->facet());
}
else
{
graph.edges[edge_index].incident_facets.insert(facet_index);
}
}
} // end for all halfedge around the facet fh
} // end if fh not yet handled
}
} // end of construct_graph_facet(...)
std::pair<bool, Halfedge_handle>
can_follow_sharp_edges(const Halfedge_handle& he) const
{
typename Geom_traits::Construct_vector_3 vector
= Geom_traits().construct_vector_3_object();
const Vertex_handle& v = he->vertex();
Halfedge_handle result = Halfedge_handle();
unsigned int nb_sharp_edges = 0;
Halfedge_around_vertex_circulator other_he = v->vertex_begin();
Halfedge_around_vertex_circulator end = other_he;
CGAL_For_all(other_he,end)
{
if(other_he->sharp())
{
++nb_sharp_edges;
if(he != other_he)
result = other_he;
}
}
CGAL_assertion(this->nb_sharp_edges(v)==nb_sharp_edges);
if(nb_sharp_edges == 2)
{
CGAL_assertion(result->vertex() == he ->vertex());
const Vector_3 vector_he = vector(he->opposite()->vertex()->point(),
he->vertex()->point());
const Vector_3 vector_result_opposite =
vector(result->vertex()->point(),
result->opposite()->vertex()->point());
return std::make_pair(true, result->opposite());
}
else
return std::make_pair(false, Halfedge_handle());
}
#ifdef CGAL_SURFACE_MESHER_POLYHEDRAL_SURFACE_USE_OCTREE
void gl_draw_facet_octree()
{
subfacets_tree_ptr->gl_draw_nodes();
}
void gl_draw_edges_octree()
{
subsegments_tree_ptr->gl_draw_nodes();
}
#endif
// draw edges
void gl_draw_sharp_edges_with_names(const float line_width,
unsigned char r,
unsigned char g,
unsigned char b)
{
::glLineWidth(line_width);
::glColor3ub(r,g,b);
for(Halfedge_iterator he = edges_begin();
he != edges_end();
he++)
{
if(he->sharp())
{
const Point_3& a = he->opposite()->vertex()->point();
const Point_3& b = he->vertex()->point();
::glPushName(he->tag()+incidence_graph.vertices.size());
::glBegin(GL_LINES);
::glVertex3d(a[0],a[1],a[2]);
::glVertex3d(b[0],b[1],b[2]);
::glEnd();
::glPopName();
}
}
}
void gl_draw_direct_triangles_with_name(bool smooth_shading,
bool use_normals,
bool inverse_normals = false)
{
// draw triangles
Facet_iterator pFacet = facets_begin();
for(;pFacet != facets_end();pFacet++)
{
::glPushName(pFacet->tag()
+ incidence_graph.vertices.size()
+ incidence_graph.edges.size());
::glBegin(GL_TRIANGLES);
gl_draw_facet(pFacet,smooth_shading,use_normals,inverse_normals);
::glEnd(); // end polygon assembly
::glPopName();
}
}
void gl_draw_almost_all_triangles(int facet_index_not_drawed,
bool smooth_shading,
bool use_normals,
bool inverse_normals = false)
{
// draw triangles
::glBegin(GL_TRIANGLES);
Facet_iterator pFacet = facets_begin();
for(;pFacet != facets_end();pFacet++)
if(pFacet->tag() != facet_index_not_drawed)
gl_draw_facet(pFacet,smooth_shading,use_normals,inverse_normals);
::glEnd(); // end polygon assembly
}
template <typename Vertex_handle>
bool vertices_not_on_same_curve(const Vertex_handle& v1,
const Vertex_handle& v2) const
{
struct Display_curves_indices {
void operator()(const std::set<int>& container) const {
std::cerr << "(";
for(std::set<int>::const_iterator
it = container.begin(),
end = container.end();
it != end;)
{
std::cerr << *it;
if(++it!=end)
std::cerr << " ";
}
std::cerr << ")";
}
};
// if(v1->point().dimension() < 0) {
// if(v2->point().dimension() < 0)
// return false; // both in volume
// else
// return true; // v1 in volume and v2 on a surface
// }
// else
// if(v2->point().dimension() < 0)
// return true; // v2 in volume, and v1 on a surface
if(v1->point().dimension() < 0 || v2->point().dimension() < 0)
return true;
std::set<int> incident_edges_v1, incident_edges_v2, intersection;
incident_edges(v1, CGAL::inserter(incident_edges_v1));
incident_edges(v2, CGAL::inserter(incident_edges_v2));
#ifdef CGAL_SURFACE_MESHER_DEBUG_INCIDES
Display_curves_indices()(incident_edges_v1);
Display_curves_indices()(incident_edges_v2);
#endif // CGAL_SURFACE_MESHER_DEBUG_INCIDES
std::set_intersection(incident_edges_v1.begin(), incident_edges_v1.end(),
incident_edges_v2.begin(), incident_edges_v2.end(),
CGAL::inserter(intersection));
return intersection.empty();
}
template <typename Vertex_handle>
bool vertices_not_on_same_surface_patch(const Vertex_handle& v1,
const Vertex_handle& v2,
const Vertex_handle& v3) const
{
// if(v1->point().dimension() < 0) {
// if(v2->point().dimension() < 0 && v3->point().dimension() < 0)
// return false; // all three vertices in volume
// else
// return true; // v1 in volume, and (v2 or v3) on a surface
// }
// else
// if(v2->point().dimension() < 0 || v3->point().dimension() < 0
if(v1->point().dimension() < 0 ||
v2->point().dimension() < 0 ||
v3->point().dimension() < 0)
return true;
std::set<int> incident_facets_v1, incident_facets_v2, incident_facets_v3, intersection;
incident_facets(v1, CGAL::inserter(incident_facets_v1));
incident_facets(v2, CGAL::inserter(incident_facets_v2));
incident_facets(v3, CGAL::inserter(incident_facets_v3));
std::set_intersection(incident_facets_v1.begin(), incident_facets_v1.end(),
incident_facets_v2.begin(), incident_facets_v2.end(),
CGAL::inserter(intersection));
if(intersection.empty())
return true;
std::set<int> intersection2;
std::set_intersection(intersection.begin(), intersection.end(),
incident_facets_v3.begin(), incident_facets_v3.end(),
CGAL::inserter(intersection2));
return intersection2.empty();
}
template <typename Vertex_handle, typename OutputIterator>
void incident_edges(const Vertex_handle& v,
OutputIterator out_it) const
{
CGAL_assertion(v->point().element_index() >= 0);
switch(v->point().dimension())
{
case 0: {
const int index = v->point().element_index();
CGAL_assertion(index < incidence_graph.vertices.size() );
CGAL_assertion(incidence_graph.vertices[index].incident_edges.size() >= 1 );
std::copy(incidence_graph.vertices[index].incident_edges.begin(),
incidence_graph.vertices[index].incident_edges.end(),
out_it);
break;
}
case 1:
*out_it++ = v->point().element_index();
break;
default:
// if dimension=2, nothing to output
// if dimension=-1 (point in volume, nothing to output)
break;
}
}
template <typename Vertex_handle, typename OutputIterator>
void incident_facets(const Vertex_handle& v,
OutputIterator out_it) const
{
CGAL_assertion(v->point().element_index() >= 0);
switch(v->point().dimension())
{
case 0: {
std::vector<int> incident_edges;
const int index = v->point().element_index();
CGAL_assertion(index < incidence_graph.vertices.size() );
std::copy(incidence_graph.vertices[index].incident_edges.begin(),
incidence_graph.vertices[index].incident_edges.end(),
std::back_inserter(incident_edges));
for(std::vector<int>::const_iterator
edge_index_it = incident_edges.begin(),
end = incident_edges.end();
edge_index_it != end;
++edge_index_it)
{
std::copy(incidence_graph.edges[*edge_index_it].incident_facets.begin(),
incidence_graph.edges[*edge_index_it].incident_facets.end(),
out_it);
}
break;
}
case 1: {
const int index = v->point().element_index();
CGAL_assertion(index < incidence_graph.edges.size() );
std::copy(incidence_graph.edges[index].incident_facets.begin(),
incidence_graph.edges[index].incident_facets.end(),
out_it);
break;
}
default: // dimension=2
*out_it++ = v->point().element_index();
break;
}
}
public:
Subfacets_tree_ptr subfacets_tree_ptr;
Subsegments_tree_ptr subsegments_tree_ptr;
typedef std::set<Vertex_handle, Compare_vertex_iterators> Input_vertices;
std::shared_ptr<Input_vertices> input_vertices_ptr;
typedef std::vector<Vertex_handle> Corner_vertices;
std::shared_ptr<Corner_vertices> corner_vertices_ptr;
typedef std::vector<Vertex_handle> Edges_vertices;
std::shared_ptr<Edges_vertices> edges_vertices_ptr;
#ifdef CGAL_SURFACE_MESHER_POLYHEDRAL_SURFACE_USE_PINPOLYHEDRON
typedef std::shared_ptr<PointInPolyhedron> PointInPolyhedron_ptr;
PointInPolyhedron_ptr pinpolyhedron_ptr;
#endif
Bbox bounding_box;
FT bounding_box_sq_radius;
double sharp_edges_angle_lower_bound;
double sharp_edges_angle_upper_bound;
double sharp_vertices_angle_lower_bound;
double sharp_vertices_angle_upper_bound;
Incidence_graph incidence_graph;
};
} // end namespace CGAL
#endif
#endif // CGAL_POLYHEDRAL_SURFACE_3_H
Reference in New Issue View Git Blame Copy Permalink