Progress on periodic optimizers

- Replaced a lot of v->point() with tr.point(v) & similar
- Fixed taking references to temporary:
  we can't do "F(construct_point_3_object()(tr.point(c, i)))" for periodic
  triangulations because tr.point(c,i) is _not_ a reference
- Added some sanity checks to debug Lloyd optimization
- Fixed some indentation
This commit is contained in:
Mael Rouxel-Labbé 2017-10-29 21:31:59 +01:00
parent 1f4475d471
commit e2b33c4948
25 changed files with 707 additions and 596 deletions

View File

@ -44,14 +44,16 @@ namespace CGAL {
//! @param c3t3 an instance of `C3T3`.
//! @param graph an instance of `TriangleMesh`.
template<class C3T3, class TriangleMesh>
void facets_in_complex_3_to_triangle_mesh(const C3T3& c3t3, TriangleMesh& graph) //complexity nlogn(number of facets on surface)
void facets_in_complex_3_to_triangle_mesh(const C3T3& c3t3, TriangleMesh& graph)
{
//complexity nlogn(number of facets on surface)
// <periodic> needs a rewrite (see medit periodic)
typedef typename boost::property_map<TriangleMesh, boost::vertex_point_t>::type VertexPointMap;
typedef typename boost::property_traits<VertexPointMap>::value_type Point_3;
typedef typename C3T3::Triangulation Tr;
typedef typename Tr::Vertex_handle Vertex_handle;
typedef typename C3T3::Triangulation Triangulation;
typedef typename Triangulation::Vertex_handle Vertex_handle;
typedef typename Triangulation::Weighted_point Weighted_point;

View File

@ -641,9 +641,6 @@ class C3T3_helpers
typedef typename Gt::Plane_3 Plane_3;
typedef typename Gt::Tetrahedron_3 Tetrahedron;
typedef typename Gt::Construct_point_3 Construct_point_3;
typedef typename Gt::Construct_weighted_point_3 Construct_weighted_point_3;
typedef typename Tr::Vertex_handle Vertex_handle;
typedef typename Tr::Facet Facet;
typedef typename Tr::Cell Cell;
@ -717,8 +714,6 @@ public:
, c3t3_(c3t3)
, tr_(c3t3.triangulation())
, domain_(domain)
, wp2p_(tr_.geom_traits().construct_point_3_object())
, p2wp_(tr_.geom_traits().construct_weighted_point_3_object())
{ }
/**
@ -1248,13 +1243,13 @@ private:
// Update facet surface center
Bare_point surface_center = CGAL::cpp0x::get<0>(intersection);
facet.first->set_facet_surface_center(facet.second,surface_center);
facet_m.first->set_facet_surface_center(facet_m.second,surface_center);
facet.first->set_facet_surface_center(facet.second, surface_center);
facet_m.first->set_facet_surface_center(facet_m.second, surface_center);
}
else
{
facet.first->set_facet_surface_center(facet.second,Bare_point());
facet_m.first->set_facet_surface_center(facet_m.second,Bare_point());
facet.first->set_facet_surface_center(facet.second, Bare_point());
facet_m.first->set_facet_surface_center(facet_m.second, Bare_point());
}
}
@ -2408,9 +2403,6 @@ private:
C3T3& c3t3_;
Tr& tr_;
const MeshDomain& domain_;
// @todo purge the two members below?
Construct_point_3 wp2p_;
Construct_weighted_point_3 p2wp_;
}; // class C3T3_helpers
@ -2425,9 +2417,10 @@ update_mesh(const Vertex_handle& old_vertex,
OutputIterator modified_vertices,
bool *could_lock_zone)
{
// std::cerr << "\nupdate_mesh[v1](" << new_position << ",\n"
// << " " << (void*)(&*old_vertex) << "=" << old_vertex->point()
// << ")\n";
// std::cerr << "\nupdate_mesh[v1](" << (void*)(&*old_vertex)
// << "=" << tr_.point(old_vertex) << ",\n"
// << " " << move << ",\n"
// << " " << new_position << ")\n";
if (could_lock_zone)
*could_lock_zone = true;
@ -2459,17 +2452,19 @@ update_mesh_no_topo_change(const Vertex_handle& old_vertex,
const Cell_vector& conflict_cells )
{
// std::cerr << "update_mesh_no_topo_change("
// << (void*)(&*old_vertex) << " = " << old_vertex->point() << ",\n"
// << (void*)(&*old_vertex) << " = " << tr_.point(old_vertex) << ",\n"
// << " " << move << ",\n"
// << " " << new_position << ")" << std::endl;
typename Gt::Construct_opposite_vector_3 cov = tr_.geom_traits().construct_opposite_vector_3_object();
//backup metadata
std::set<Cell_data_backup> cells_backup;
fill_cells_backup(conflict_cells, cells_backup);
// Get old values
criterion.before_move(c3t3_cells(conflict_cells));
Weighted_point old_position = old_vertex->point();
Weighted_point old_position = tr_.point(old_vertex); // intentional copy
// Move point
reset_circumcenter_cache(conflict_cells);
@ -2496,10 +2491,6 @@ update_mesh_no_topo_change(const Vertex_handle& old_vertex,
//sliver caches have been updated by valid_move
reset_sliver_cache(conflict_cells);
typename Gt::Construct_opposite_vector_3 cov =
tr_.geom_traits().construct_opposite_vector_3_object();
move_point_no_topo_change(old_vertex, cov(move), old_position);
//restore meta-data (cells should have same connectivity as before move)
@ -2507,7 +2498,7 @@ update_mesh_no_topo_change(const Vertex_handle& old_vertex,
CGAL_assertion(conflict_cells.size() >= cells_backup.size());
restore_from_cells_backup(conflict_cells, cells_backup);
return std::make_pair(false,old_vertex);
return std::make_pair(false, old_vertex);
}
}
@ -2522,8 +2513,8 @@ update_mesh_topo_change(const Vertex_handle& old_vertex,
OutputIterator modified_vertices,
bool *could_lock_zone)
{
// std::cerr << "update_mesh_topo_change("
// << (void*)(&*old_vertex) << "=" << old_vertex->point()
// std::cerr << "update_mesh_topo_change(" << (void*)(&*old_vertex)
// << "=" << tr_.point(old_vertex)
// << " " << move << ",\n"
// << " " << new_position << ",\n"
// << ")" << std::endl;
@ -2558,7 +2549,7 @@ update_mesh_topo_change(const Vertex_handle& old_vertex,
CGAL_assertion(conflict_cells.size() == cells_backup.size());
criterion.before_move(c3t3_cells(conflict_cells));
Weighted_point old_position = old_vertex->point();
Weighted_point old_position = tr_.point(old_vertex); // intentional copy
// Keep old boundary
Vertex_set old_incident_surface_vertices;
@ -2581,7 +2572,7 @@ update_mesh_topo_change(const Vertex_handle& old_vertex,
CGAL::Emptyset_iterator());
// If nothing changed, return
if ( old_position == new_vertex->point() )
if ( old_position == tr_.point(new_vertex) )
{
// std::cerr << "update_mesh_topo_change: no move!\n";
// check_c3t3(c3t3_);
@ -2639,9 +2630,9 @@ update_mesh(const Vertex_handle& old_vertex,
OutputIterator modified_vertices,
bool fill_vertices)
{
// std::cerr << "\nupdate_mesh[v2](" << new_position << ",\n"
// << " " << (void*)(&*old_vertex) << "=" << old_vertex->point()
// << ")\n";
// std::cerr << "\nupdate_mesh[v2](" << (void*)(&*old_vertex)
// << "=" << tr_.point(old_vertex) << ",\n"
// << " " << new_position << ")\n";
Cell_vector outdated_cells;
Vertex_handle new_vertex = move_point(old_vertex, move,
@ -2670,7 +2661,9 @@ C3T3_helpers<C3T3,MD>::
rebuild_restricted_delaunay(OutdatedCells& outdated_cells,
Moving_vertices_set& moving_vertices)
{
typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object();
typename Gt::Equal_3 equal = tr_.geom_traits().equal_3_object();
typename OutdatedCells::iterator first_cell = outdated_cells.begin();
typename OutdatedCells::iterator last_cell = outdated_cells.end();
Update_c3t3 updater(domain_,c3t3_);
@ -2768,8 +2761,8 @@ rebuild_restricted_delaunay(OutdatedCells& outdated_cells,
it != vertex_to_proj.end() ;
++it )
{
Bare_point initial_position = wp2p_((*it)->point());
Bare_point new_pos = project_on_surface(*it, initial_position);
const Weighted_point& initial_position = tr_.point(*it);
Bare_point new_pos = project_on_surface(*it, cp(initial_position));
std::cout << "project: " << initial_position << " got: " << new_pos << std::endl;
@ -2779,7 +2772,7 @@ rebuild_restricted_delaunay(OutdatedCells& outdated_cells,
// Update moving vertices (it becomes new_vertex)
moving_vertices.erase(*it);
Vector_3 move(initial_position, new_pos);
Vector_3 move(cp(initial_position), new_pos);
std::cout << "projection move: " << move << std::endl;
CGAL_assertion(CGAL::abs(move.x()) < 0.8*(tr_.domain().xmax() - tr_.domain().xmin()));
@ -2807,8 +2800,8 @@ rebuild_restricted_delaunay(ForwardIterator first_cell,
ForwardIterator last_cell,
Moving_vertices_set& moving_vertices)
{
typename Gt::Equal_3 equal = tr_.geom_traits().equal_3_object();
typename Gt::Construct_vector_3 vector = tr_.geom_traits().construct_vector_3_object();
typename Gt::Equal_3 equal = tr_.geom_traits().equal_3_object();
Update_c3t3 updater(domain_,c3t3_);
@ -2886,16 +2879,17 @@ rebuild_restricted_delaunay(ForwardIterator first_cell,
end = vertex_to_proj.end();
for ( ; it != end; ++it )
{
const Bare_point& initial_pos = wp2p_((it->first)->point());
Bare_point new_pos = project_on_surface(it->first, initial_pos, it->second);
Vertex_handle vh = it->first;
const Weighted_point& initial_position = tr_.point(vh);
Bare_point new_pos = project_on_surface(vh, cp(initial_position), it->second);
if ( ! equal(new_pos, Bare_point()) )
{
//freezing needs 'erase' to be done before the vertex is actually destroyed
// Update moving vertices (it becomes new_vertex)
moving_vertices.erase(it->first);
moving_vertices.erase(vh);
Vertex_handle new_vertex = update_mesh(it->first, vector(initial_pos, new_pos));
Vertex_handle new_vertex = update_mesh(vh, vector(cp(initial_position), new_pos));
c3t3_.set_dimension(new_vertex, 2);
moving_vertices.insert(new_vertex);
@ -2925,22 +2919,20 @@ move_point(const Vertex_handle& old_vertex,
OutdatedCellsOutputIterator outdated_cells,
DeletedCellsOutputIterator deleted_cells) const
{
// std::cerr << "C3T3_helpers::move_point[v2]("
// << (void*)(&*old_vertex) << " = " << old_vertex->point()
// << " , move: " << move << ")\n";
// std::cerr << "C3T3_helpers::move_point[v2](" << (void*)(&*old_vertex)
// << " = " << tr_.point(old_vertex) << ",\n"
// << " " << move << ")\n";
typename Gt::Construct_translated_point_3 translate =
tr_.geom_traits().construct_translated_point_3_object();
typename Gt::Construct_weighted_point_3 p2wp =
tr_.geom_traits().construct_weighted_point_3_object();
typename Gt::Construct_point_3 wp2p =
tr_.geom_traits().construct_point_3_object();
typename Gt::Construct_translated_point_3 translate = tr_.geom_traits().construct_translated_point_3_object();
typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object();
typename Gt::Construct_weighted_point_3 cwp = tr_.geom_traits().construct_weighted_point_3_object();
Cell_vector incident_cells_;
incident_cells_.reserve(64);
tr_.incident_cells(old_vertex, std::back_inserter(incident_cells_));
const Weighted_point& new_position = p2wp(translate(wp2p(tr_.point(old_vertex)), move));
const Weighted_point& position = tr_.point(old_vertex);
const Weighted_point& new_position = cwp(translate(cp(position), move));
if ( Th().no_topological_change(tr_, old_vertex, move, new_position, incident_cells_) )
{
@ -2964,18 +2956,16 @@ move_point(const Vertex_handle& old_vertex,
Outdated_cell_set& outdated_cells_set,
Moving_vertices_set& moving_vertices) const
{
typename Gt::Construct_translated_point_3 translate =
tr_.geom_traits().construct_translated_point_3_object();
typename Gt::Construct_weighted_point_3 p2wp =
tr_.geom_traits().construct_weighted_point_3_object();
typename Gt::Construct_point_3 wp2p =
tr_.geom_traits().construct_point_3_object();
typename Gt::Construct_translated_point_3 translate = tr_.geom_traits().construct_translated_point_3_object();
typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object();
typename Gt::Construct_weighted_point_3 cwp = tr_.geom_traits().construct_weighted_point_3_object();
Cell_vector incident_cells_;
incident_cells_.reserve(64);
tr_.incident_cells(old_vertex, std::back_inserter(incident_cells_));
const Weighted_point& new_position = p2wp(translate(wp2p(tr_.point(old_vertex)), move));
const Weighted_point& position = tr_.point(old_vertex);
const Weighted_point& new_position = cwp(translate(cp(position), move));
if ( Th().no_topological_change(tr_, old_vertex, move, new_position, incident_cells_) )
{
@ -3029,14 +3019,12 @@ move_point(const Vertex_handle& old_vertex,
}
//======= /Get incident cells ==========
typename Gt::Construct_translated_point_3 translate =
tr_.geom_traits().construct_translated_point_3_object();
typename Gt::Construct_weighted_point_3 p2wp =
tr_.geom_traits().construct_weighted_point_3_object();
typename Gt::Construct_point_3 wp2p =
tr_.geom_traits().construct_point_3_object();
typename Gt::Construct_translated_point_3 translate = tr_.geom_traits().construct_translated_point_3_object();
typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object();
typename Gt::Construct_weighted_point_3 cwp = tr_.geom_traits().construct_weighted_point_3_object();
const Weighted_point& new_position = p2wp(translate(wp2p(tr_.point(old_vertex)), move));
const Weighted_point& position = tr_.point(old_vertex);
const Weighted_point& new_position = cwp(translate(cp(position), move));
if (!try_lock_point(new_position)) // LOCK
{
@ -3198,7 +3186,7 @@ move_point_topo_change_conflict_zone_known(
//o.w. deleted_cells will point to null pointer or so and crash
const
{
Weighted_point old_position = old_vertex->point();
Weighted_point old_position = tr_.point(old_vertex); // intentional copy
// make one set with conflict zone
Cell_set conflict_zone;
std::set_union(insertion_conflict_cells_begin, insertion_conflict_cells_end,
@ -3325,16 +3313,11 @@ project_on_surface_aux(const Bare_point& p,
typedef typename Gt::Segment_3 Segment_3;
// Build a segment directed as projection_direction,
typename Gt::Compute_squared_distance_3 sq_distance =
tr_.geom_traits().compute_squared_distance_3_object();
typename Gt::Compute_squared_length_3 sq_length =
tr_.geom_traits().compute_squared_length_3_object();
typename Gt::Construct_scaled_vector_3 scale =
tr_.geom_traits().construct_scaled_vector_3_object();
typename Gt::Is_degenerate_3 is_degenerate =
tr_.geom_traits().is_degenerate_3_object();
typename Gt::Construct_translated_point_3 translate =
tr_.geom_traits().construct_translated_point_3_object();
typename Gt::Compute_squared_distance_3 sq_distance = tr_.geom_traits().compute_squared_distance_3_object();
typename Gt::Compute_squared_length_3 sq_length = tr_.geom_traits().compute_squared_length_3_object();
typename Gt::Construct_scaled_vector_3 scale = tr_.geom_traits().construct_scaled_vector_3_object();
typename Gt::Is_degenerate_3 is_degenerate = tr_.geom_traits().is_degenerate_3_object();
typename Gt::Construct_translated_point_3 translate = tr_.geom_traits().construct_translated_point_3_object();
typename MD::Construct_intersection construct_intersection =
domain_.construct_intersection_object();
@ -3388,6 +3371,8 @@ get_least_square_surface_plane(const Vertex_handle& v,
Bare_point& reference_point,
Surface_patch_index patch_index) const
{
typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object();
// Get incident facets
Facet_vector facets;
# ifdef CGAL_LINKED_WITH_TBB
@ -3416,7 +3401,8 @@ get_least_square_surface_plane(const Vertex_handle& v,
const int& i = fit->second;
// @fixme really ugly
const Bare_point& bp = tr_.get_closest_point(wp2p_(tr_.point(v)),
const Weighted_point& position = tr_.point(v);
const Bare_point& bp = tr_.get_closest_point(cp(position),
cell->get_facet_surface_center(i));
surface_point_vector.push_back(bp);
}
@ -3453,6 +3439,7 @@ project_on_surface(const Vertex_handle& v,
{
// return domain_.project_on_surface(p);
typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object();
typename Gt::Equal_3 equal = tr_.geom_traits().equal_3_object();
// Get plane
@ -3463,8 +3450,9 @@ project_on_surface(const Vertex_handle& v,
return p;
// Project
if ( ! equal(p, wp2p_(v->point())) )
return project_on_surface_aux(p, wp2p_(v->point()), plane.orthogonal_vector());
const Weighted_point& position = tr_.point(v);
if ( ! equal(p, cp(position)) )
return project_on_surface_aux(p, cp(position), plane.orthogonal_vector());
else
return project_on_surface_aux(p, reference_point, plane.orthogonal_vector());
}

View File

@ -66,7 +66,7 @@ public:
: Base(ch)
{
const Weighted_point& p = ch->vertex(0)->point();
const Weighted_point& p = ch->vertex(0)->point(); // <periodic> ?
const Weighted_point& q = ch->vertex(1)->point();
const Weighted_point& r = ch->vertex(2)->point();
const Weighted_point& s = ch->vertex(3)->point();

View File

@ -44,7 +44,8 @@ template <typename C3t3,
Output_rep<typename C3t3::Subdomain_index>::is_specialized)
>
struct Dump_c3t3 {
void dump_c3t3(const C3t3& c3t3, std::string prefix) const {
void dump_c3t3(const C3t3& c3t3, std::string prefix) const
{
std::clog<<"======dump c3t3===== to: " << prefix << std::endl;
std::ofstream medit_file((prefix+".mesh").c_str());
medit_file.precision(17);
@ -64,7 +65,8 @@ struct Dump_c3t3 {
}; // end struct template Dump_c3t3<C3t3, bool>
template <typename C3t3>
struct Dump_c3t3<C3t3, false> {
struct Dump_c3t3<C3t3, false>
{
void dump_c3t3(const C3t3&, std::string) {
std::cerr << "Warning " << __FILE__ << ":" << __LINE__ << "\n"
<< " the c3t3 object of following type:\n"
@ -104,8 +106,9 @@ struct Dump_c3t3<C3t3, false> {
}; // end struct template specialization Dump_c3t3<C3t3, false>
template <typename C3t3>
void dump_c3t3_edges(const C3t3& c3t3, std::string prefix) {
typename C3t3::Triangulation::Geom_traits::Construct_point_3 wp2p =
void dump_c3t3_edges(const C3t3& c3t3, std::string prefix)
{
typename C3t3::Triangulation::Geom_traits::Construct_point_3 cp =
c3t3.triangulation().geom_traits().construct_point_3_object();
std::ofstream file((prefix+".polylines.txt").c_str());
@ -118,12 +121,14 @@ void dump_c3t3_edges(const C3t3& c3t3, std::string prefix) {
const typename C3t3::Triangulation::Cell_handle c = edge_it->first;
const int i = edge_it->second;
const int j = edge_it->third;
file << "2 " << wp2p(c->vertex(i)->point())
<< " " << wp2p(c->vertex(j)->point()) << "\n";
const typename C3t3::Weighted_point& ei = c3t3.triangulation().point(c, i);
const typename C3t3::Weighted_point& ej = c3t3.triangulation().point(c, j);
file << "2 " << cp(ei) << " " << cp(ej) << "\n";
}
}
template <typename C3t3>
void dump_c3t3(const C3t3& c3t3, std::string prefix) {
void dump_c3t3(const C3t3& c3t3, std::string prefix)
{
if(!prefix.empty()) {
Dump_c3t3<C3t3> dump;
dump.dump_c3t3(c3t3, prefix);

View File

@ -78,7 +78,7 @@ public:
Squared_radius_orthogonal_sphere sq_radius_ortho_sphere =
Gt().compute_squared_radius_smallest_orthogonal_sphere_3_object();
Weighted_point p1 = fh.first->vertex ((fh.second+1)&3)->point();
Weighted_point p1 = fh.first->vertex ((fh.second+1)&3)->point(); // <periodic> ?
Weighted_point p2 = fh.first->vertex ((fh.second+2)&3)->point();
Weighted_point p3 = fh.first->vertex ((fh.second+3)&3)->point();

View File

@ -55,6 +55,7 @@ class Lloyd_move
typedef typename Tr::Cell_handle Cell_handle;
typedef typename Tr::Bare_point Bare_point;
typedef typename Tr::Weighted_point Weighted_point;
typedef typename std::vector<Facet> Facet_vector;
typedef typename std::vector<Cell_handle> Cell_vector;
@ -79,7 +80,7 @@ public:
const Sizing_field& sizing_field = Sizing_field() ) const
{
// std::cout << "computing move of: " << &*v
// << " pos: " << v->point().point()
// << " pos: " << c3t3.triangulation().point(v)
// << " dim: " << c3t3.in_dimension(v) << std::endl;
switch ( c3t3.in_dimension(v) )
@ -260,7 +261,7 @@ private:
{
const Tr& tr = c3t3.triangulation();
typename Gt::Construct_point_3 wp2p = tr.geom_traits().construct_point_3_object();
typename Gt::Construct_point_3 cp = tr.geom_traits().construct_point_3_object();
Facet_vector incident_facets;
incident_facets.reserve(64);
@ -281,9 +282,10 @@ private:
const Cell_handle& cell = fit->first;
const int& i = fit->second;
// @fixme really ugly
points.push_back(tr.get_closest_point(wp2p(tr.point(v)),
cell->get_facet_surface_center(i)));
// @fixme 'get_closest_point' is really ugly
const Weighted_point& position = tr.point(v);
points.push_back(tr.get_closest_point(cp(position),
cell->get_facet_surface_center(i)));
}
return points;
@ -298,20 +300,19 @@ private:
const C3T3& c3t3,
const Sizing_field& sizing_field) const
{
typename Gt::Construct_point_3 cp =
c3t3.triangulation().geom_traits().construct_point_3_object();
typename Gt::Construct_vector_3 vector =
c3t3.triangulation().geom_traits().construct_vector_3_object();
typename Gt::Construct_point_3 cp = c3t3.triangulation().geom_traits().construct_point_3_object();
typename Gt::Construct_vector_3 vector = c3t3.triangulation().geom_traits().construct_vector_3_object();
const Bare_point& p = cp(v->point());
const Weighted_point position = c3t3.triangulation().point(v);
const Bare_point& p = cp(position);
FT da = density_1d(a,v,sizing_field);
FT db = density_1d(b,v,sizing_field);
FT da = density_1d(a, v, sizing_field);
FT db = density_1d(b, v, sizing_field);
CGAL_assertion( !is_zero(da + db) );
return ( (vector(p,a)*da + vector(p,b)*db) / (da + db) );
}
/**
* Return the move from \c v to the centroid of triangle [a,b,c].
*/
@ -322,12 +323,11 @@ private:
const C3T3& c3t3,
const Sizing_field& sizing_field) const
{
typename Gt::Construct_vector_3 vector =
c3t3.triangulation().geom_traits().construct_vector_3_object();
typename Gt::Construct_point_3 cp =
c3t3.triangulation().geom_traits().construct_point_3_object();
typename Gt::Construct_point_3 cp = c3t3.triangulation().geom_traits().construct_point_3_object();
typename Gt::Construct_vector_3 vector = c3t3.triangulation().geom_traits().construct_vector_3_object();
const Bare_point& p = cp(v->point());
const Weighted_point& position = c3t3.triangulation().point(v);
const Bare_point& p = cp(position);
FT da = density_2d<true>(a,v,sizing_field);
FT db = density_2d<false>(b,v,sizing_field);
@ -360,7 +360,7 @@ private:
// Project all points to the plane
std::transform(first, last, first, Project_on_plane(plane, c3t3));
CGAL_assertion(std::distance(first,last) > 3);
CGAL_assertion(std::distance(first, last) > 3);
// Get 2D-3D transformations
Aff_transformation_3 to_3d = compute_to_3d_transform(plane, *first, c3t3);
@ -403,17 +403,14 @@ private:
{
CGAL_precondition(std::distance(first,last) >= 3);
typename Gt::Construct_vector_3 vector =
c3t3.triangulation().geom_traits().construct_vector_3_object();
typename Gt::Construct_centroid_3 centroid =
c3t3.triangulation().geom_traits().construct_centroid_3_object();
typename Gt::Compute_area_3 area =
c3t3.triangulation().geom_traits().compute_area_3_object();
typename Gt::Construct_point_3 cp =
c3t3.triangulation().geom_traits().construct_point_3_object();
typename Gt::Compute_area_3 area = c3t3.triangulation().geom_traits().compute_area_3_object();
typename Gt::Construct_centroid_3 centroid = c3t3.triangulation().geom_traits().construct_centroid_3_object();
typename Gt::Construct_point_3 cp = c3t3.triangulation().geom_traits().construct_point_3_object();
typename Gt::Construct_vector_3 vector = c3t3.triangulation().geom_traits().construct_vector_3_object();
// Vertex current position
const Bare_point& vertex_position = cp(v->point());
const Weighted_point& vertex_weighted_position = c3t3.triangulation().point(v);
const Bare_point& vertex_position = cp(vertex_weighted_position);
// Use as reference point to triangulate
const Bare_point& a = *first++;
@ -457,10 +454,8 @@ private:
const Bare_point& reference_point,
const C3T3& c3t3) const
{
typename Gt::Construct_base_vector_3 base =
c3t3.triangulation().geom_traits().construct_base_vector_3_object();
typename Gt::Construct_orthogonal_vector_3 orthogonal_vector =
c3t3.triangulation().geom_traits().construct_orthogonal_vector_3_object();
typename Gt::Construct_base_vector_3 base = c3t3.triangulation().geom_traits().construct_base_vector_3_object();
typename Gt::Construct_orthogonal_vector_3 orthogonal_vector = c3t3.triangulation().geom_traits().construct_orthogonal_vector_3_object();
Vector_3 u = base(plane, 1);
u = u / CGAL::sqrt(u*u);
@ -553,19 +548,20 @@ private:
CGAL_assertion(c3t3.is_in_complex(current_cell));
// a & b are fixed points
const Bare_point& a = cp(v->point());
const Weighted_point& wa = tr.point(v);
const Bare_point& a = cp(wa);
Bare_point b = tr.dual(current_cell);
const Bare_point& a_b = cp(tr.point(current_cell, current_cell->index(v)));
Vector_3 ba = Vector_3(a_b, b);
const Weighted_point& a_b = tr.point(current_cell, current_cell->index(v));
Vector_3 ba = Vector_3(cp(a_b), b);
++current_cell;
CGAL_assertion(c3t3.is_in_complex(current_cell));
CGAL_assertion(current_cell != done);
// c & d are moving points
Bare_point c = tr.dual(current_cell);
const Bare_point& a_c = cp(tr.point(current_cell, current_cell->index(v)));
Vector_3 ca = Vector_3(a_c, c);
const Weighted_point& a_c = tr.point(current_cell, current_cell->index(v));
Vector_3 ca = Vector_3(cp(a_c), c);
++current_cell;
CGAL_assertion(current_cell != done);
@ -573,8 +569,8 @@ private:
{
CGAL_assertion(c3t3.is_in_complex(current_cell));
Bare_point d = tr.dual(current_cell);
const Bare_point& a_d = cp(tr.point(current_cell, current_cell->index(v)));
Vector_3 da = Vector_3(a_d, d);
const Weighted_point& a_d = tr.point(current_cell, current_cell->index(v));
Vector_3 da = Vector_3(cp(a_d), d);
Tetrahedron_3 tet = tetrahedron(a, translate(a, ba), translate(a, ca), translate(a, da));
Bare_point tet_centroid = centroid(tet);

View File

@ -291,7 +291,7 @@ public:
default :
#ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS
std::cerr << "singular edge...\n";
std::cerr << v->point() << std::endl;
std::cerr << tr_.point(v) << std::endl;
#endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS
return SINGULAR;
}
@ -302,13 +302,13 @@ public:
if(nb_components > 1) {
#ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS
std::cerr << "singular vertex: nb_components=" << nb_components << std::endl;
std::cerr << v->point() << std::endl;
std::cerr << tr_.point(v) << std::endl;
#endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS
return SINGULAR;
}
else { // REGULAR OR BOUNDARY
#ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS
std::cerr << "regular or boundary: " << v->point() << std::endl;
std::cerr << "regular or boundary: " << tr_.point(v) << std::endl;
#endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS
if (number_of_boundary_incident_edges != 0)
return BOUNDARY;
@ -933,7 +933,7 @@ Mesh_complex_3_in_triangulation_3_base<Tr,Ct>::add_to_complex(
if (j != i) {
#ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS
if(cell->vertex(j)->is_c2t3_cache_valid())
std::cerr << "(" << cell->vertex(j)->point() << ")->invalidate_c2t3_cache()\n";
std::cerr << "(" << tr_.point(cell, j) << ")->invalidate_c2t3_cache()\n";
#endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS
cell->vertex(j)->invalidate_c2t3_cache();
}
@ -984,7 +984,7 @@ Mesh_complex_3_in_triangulation_3_base<Tr,Ct>::remove_from_complex(const Facet&
if (j != facet.second) {
#ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS
if(cell->vertex(j)->is_c2t3_cache_valid())
std::cerr << "(" << cell->vertex(j)->point() << ")->invalidate_c2t3_cache()\n";
std::cerr << "(" << tr_.point(cell, j) << ")->invalidate_c2t3_cache()\n";
#endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS
cell->vertex(j)->invalidate_c2t3_cache();
}
@ -1008,12 +1008,12 @@ bbox() const
}
typename Tr::Finite_vertices_iterator vit = tr_.finite_vertices_begin();
Bbox_3 result = (vit++)->point().bbox();
Bbox_3 result = tr_.point(vit++).bbox();
for(typename Tr::Finite_vertices_iterator end = tr_.finite_vertices_end();
vit != end ; ++vit)
{
result = result + vit->point().bbox();
result = result + tr_.point(vit).bbox();
}
return result;

View File

@ -1,4 +1,4 @@
// Copyright (c) 2009 INRIA Sophia-Antipolis (France).
// Copyright (c) 2009 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
@ -398,12 +398,9 @@ private:
// operator()
void operator()(const Vertex_handle& oldv) const
{
typename Gt::Construct_point_3 wp2p =
m_gt.construct_point_3_object();
typename Gt::Construct_weighted_point_3 p2wp =
m_gt.construct_weighted_point_3_object();
typename Gt::Construct_translated_point_3 translate =
m_gt.construct_translated_point_3_object();
typename Gt::Construct_point_3 cp = m_gt.construct_point_3_object();
typename Gt::Construct_weighted_point_3 cwp = m_gt.construct_weighted_point_3_object();
typename Gt::Construct_translated_point_3 translate = m_gt.construct_translated_point_3_object();
Vector_3 move = m_mgo.compute_move(oldv);
if ( CGAL::NULL_VECTOR != move )
@ -412,7 +409,8 @@ private:
if(MGO::Sizing_field::is_vertex_update_needed)
{
Bare_point new_position = translate(wp2p(oldv->point()), move);
const Weighted_point& position = m_gt.tr_.point(oldv);
Bare_point new_position = translate(cp(position), move);
size = m_sizing_field(new_position, oldv);
}
@ -464,13 +462,13 @@ private:
// operator()
void operator()(Vertex& v) const
{
typename Gt::Construct_point_3 wp2p =
m_gt.construct_point_3_object();
typename Gt::Construct_point_3 cp = m_gt.construct_point_3_object();
Vertex_handle vh
= Tr::Triangulation_data_structure::Vertex_range::s_iterator_to(v);
const Weighted_point& position = m_mgo.tr_.point(vh);
m_local_lists.local().push_back(
std::make_pair(wp2p(v.point()), m_mgo.average_circumradius_length(vh)));
std::make_pair(cp(position), m_mgo.average_circumradius_length(vh)));
}
};
@ -508,7 +506,7 @@ private:
void operator()( const tbb::blocked_range<size_t>& r ) const
{
typename Gt::Construct_translated_point_3 translate =
m_gt.construct_translated_point_3_object();
m_mgo.tr_.geom_traits().construct_translated_point_3_object();
for( size_t i = r.begin() ; i != r.end() ; ++i)
{
@ -816,10 +814,8 @@ compute_moves(Moving_vertices_set& moving_vertices)
else
#endif // CGAL_LINKED_WITH_TBB
{
typename Gt::Construct_point_3 wp2p =
tr_.geom_traits().construct_point_3_object();
typename Gt::Construct_translated_point_3 translate =
tr_.geom_traits().construct_translated_point_3_object();
typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object();
typename Gt::Construct_translated_point_3 translate = tr_.geom_traits().construct_translated_point_3_object();
// Get move for each moving vertex
typename Moving_vertices_set::iterator vit = moving_vertices.begin();
@ -834,7 +830,8 @@ compute_moves(Moving_vertices_set& moving_vertices)
FT size = 0.;
if(Sizing_field::is_vertex_update_needed)
{
Bare_point new_position = translate(wp2p(tr_.point(oldv)), move);
const Weighted_point& position = tr_.point(oldv);
Bare_point new_position = translate(cp(position), move);
// std::cout << "new position: " << new_position << std::endl;
size = sizing_field_(new_position, oldv);
}
@ -870,14 +867,10 @@ typename Mesh_global_optimizer<C3T3,Md,Mf,V_>::Vector_3
Mesh_global_optimizer<C3T3,Md,Mf,V_>::
compute_move(const Vertex_handle& v)
{
typename Gt::Compute_squared_length_3 sq_length =
tr_.geom_traits().compute_squared_length_3_object();
typename Gt::Construct_vector_3 vector =
tr_.geom_traits().construct_vector_3_object();
typename Gt::Construct_translated_point_3 translate =
tr_.geom_traits().construct_translated_point_3_object();
typename Gt::Construct_point_3 wp2p =
tr_.geom_traits().construct_point_3_object();
typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object();
typename Gt::Compute_squared_length_3 sq_length = tr_.geom_traits().compute_squared_length_3_object();
typename Gt::Construct_translated_point_3 translate = tr_.geom_traits().construct_translated_point_3_object();
typename Gt::Construct_vector_3 vector = tr_.geom_traits().construct_vector_3_object();
Cell_vector incident_cells;
incident_cells.reserve(64);
@ -904,16 +897,17 @@ compute_move(const Vertex_handle& v)
// Project surface vertex
if ( c3t3_.in_dimension(v) == 2 )
{
Bare_point new_position = translate(wp2p(tr_.point(v)), move);
const Weighted_point& position = tr_.point(v);
Bare_point new_position = translate(cp(position), move);
// std::cout << "new pos before projection: " << new_position << std::endl;
Bare_point projected_new_position = helper_.project_on_surface(v, new_position);
// std::cout << "projected: " << projected_new_position << std::endl;
move = vector(wp2p(tr_.point(v)), projected_new_position);
move = vector(cp(position), projected_new_position);
}
FT local_move_sq_ratio = sq_length(move) / local_sq_size;
// std::cout << "optimizing point: " << tr_.point(v) << std::endl;
// std::cout << "optimizing point: " << position << std::endl;
// std::cout << "moving with: " << move << std::endl;
CGAL_assertion(CGAL::abs(move.x()) < 0.8*(tr_.domain().xmax() - tr_.domain().xmin()));
@ -1074,16 +1068,15 @@ fill_sizing_field()
else
#endif //CGAL_LINKED_WITH_TBB
{
typename Gt::Construct_point_3 wp2p =
tr_.geom_traits().construct_point_3_object();
typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object();
// Fill map with local size
for(typename Tr::Finite_vertices_iterator vit = tr_.finite_vertices_begin();
vit != tr_.finite_vertices_end();
++vit)
{
value_map.insert(std::make_pair(wp2p(tr_.point(vit)),
average_circumradius_length(vit)));
const Weighted_point& position = tr_.point(vit);
value_map.insert(std::make_pair(cp(position), average_circumradius_length(vit)));
}
}
@ -1209,14 +1202,13 @@ typename Mesh_global_optimizer<C3T3,Md,Mf,V_>::FT
Mesh_global_optimizer<C3T3,Md,Mf,V_>::
sq_circumradius_length(const Cell_handle& cell, const Vertex_handle& v) const
{
typename Gt::Compute_squared_distance_3 sq_distance =
tr_.geom_traits().compute_squared_distance_3_object();
typename Gt::Construct_point_3 wp2p =
tr_.geom_traits().construct_point_3_object();
typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object();
typename Gt::Compute_squared_distance_3 sq_distance = tr_.geom_traits().compute_squared_distance_3_object();
const Bare_point circumcenter = tr_.dual(cell);
const Weighted_point& position = tr_.point(cell, cell->index(v));
return ( sq_distance(wp2p(tr_.point(cell, cell->index(v))), circumcenter) );
return ( sq_distance(cp(position), circumcenter) );
}
} // end namespace Mesh_3

View File

@ -95,12 +95,13 @@ class Mesh_sizing_field
typename Tr::Concurrency_tag>
{
// Types
typedef typename Tr::Geom_traits Gt;
typedef typename Tr::Bare_point Bare_point;
typedef typename Gt::FT FT;
typedef typename Tr::Geom_traits Gt;
typedef typename Tr::Bare_point Bare_point;
typedef typename Tr::Weighted_point Weighted_point;
typedef typename Gt::FT FT;
typedef typename Tr::Vertex_handle Vertex_handle;
typedef typename Tr::Cell_handle Cell_handle;
typedef typename Tr::Vertex_handle Vertex_handle;
typedef typename Tr::Cell_handle Cell_handle;
public:
// update vertices of mesh triangulation ?
@ -175,12 +176,14 @@ fill(const std::map<Bare_point, FT>& value_map)
{
typedef typename Tr::Finite_vertices_iterator Fvi;
for ( Fvi vit = tr_.finite_vertices_begin() ;
vit != tr_.finite_vertices_end() ;
++ vit )
typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object();
for ( Fvi vit = tr_.finite_vertices_begin(); vit != tr_.finite_vertices_end(); ++ vit )
{
const Weighted_point& position = tr_.point(vit);
typename std::map<Bare_point, FT>::const_iterator find_result =
value_map.find(tr_.geom_traits().construct_point_3_object()(vit->point()));
value_map.find(cp(position));
if ( find_result != value_map.end() )
{
@ -200,17 +203,16 @@ typename Mesh_sizing_field<Tr,B>::FT
Mesh_sizing_field<Tr,B>::
operator()(const Bare_point& p, const Cell_handle& c) const
{
typename Gt::Construct_weighted_point_3 p2wp =
tr_.geom_traits().construct_weighted_point_3_object();
typename Gt::Construct_weighted_point_3 cwp = tr_.geom_traits().construct_weighted_point_3_object();
#ifdef CGAL_MESH_3_SIZING_FIELD_INEXACT_LOCATE
//use the inexact locate (much faster than locate) to get a hint
//and then use locate to check whether p is really inside hint
// if not, an exact locate will be performed
Cell_handle hint = tr_.inexact_locate(p2wp(p),c);
const Cell_handle cell = tr_.locate(p2wp(p), hint);
Cell_handle hint = tr_.inexact_locate(cwp(p),c);
const Cell_handle cell = tr_.locate(cwp(p), hint);
#else
const Cell_handle cell = tr_.locate(p2wp(p),c);
const Cell_handle cell = tr_.locate(cwp(p),c);
#endif
this->set_last_cell(cell);
@ -244,8 +246,8 @@ typename Mesh_sizing_field<Tr,B>::FT
Mesh_sizing_field<Tr,B>::
interpolate_on_cell_vertices(const Bare_point& p, const Cell_handle& cell) const
{
typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object();
typename Gt::Compute_volume_3 volume = tr_.geom_traits().compute_volume_3_object();
typename Gt::Construct_point_3 wp2p = tr_.geom_traits().construct_point_3_object();
// Interpolate value using tet vertices values
const FT& va = cell->vertex(0)->meshing_info();
@ -253,10 +255,14 @@ interpolate_on_cell_vertices(const Bare_point& p, const Cell_handle& cell) const
const FT& vc = cell->vertex(2)->meshing_info();
const FT& vd = cell->vertex(3)->meshing_info();
const Bare_point& a = wp2p(tr_.point(cell, 0));
const Bare_point& b = wp2p(tr_.point(cell, 1));
const Bare_point& c = wp2p(tr_.point(cell, 2));
const Bare_point& d = wp2p(tr_.point(cell, 3));
const Weighted_point& wa = tr_.point(cell, 0);
const Weighted_point& wb = tr_.point(cell, 1);
const Weighted_point& wc = tr_.point(cell, 2);
const Weighted_point& wd = tr_.point(cell, 3);
const Bare_point& a = cp(wa);
const Bare_point& b = cp(wb);
const Bare_point& c = cp(wc);
const Bare_point& d = cp(wd);
const FT abcp = CGAL::abs(volume(a,b,c,p));
const FT abdp = CGAL::abs(volume(a,d,b,p));
@ -265,13 +271,12 @@ interpolate_on_cell_vertices(const Bare_point& p, const Cell_handle& cell) const
// If volume is 0, then compute the average value
if ( is_zero(abcp+abdp+acdp+bcdp) )
return (va+vb+vc+vd)/4.;
return (va+vb+vc+vd) / 4.;
return ( (abcp*vd + abdp*vc + acdp*vb + bcdp*va) / (abcp+abdp+acdp+bcdp) );
}
template <typename Tr, bool B>
typename Mesh_sizing_field<Tr,B>::FT
Mesh_sizing_field<Tr,B>::
@ -279,7 +284,7 @@ interpolate_on_facet_vertices(const Bare_point& p, const Cell_handle& cell) cons
{
typename Gt::Compute_area_3 area = tr_.geom_traits().compute_area_3_object();
typename Gt::Construct_point_3 wp2p = tr_.geom_traits().construct_point_3_object();
typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object();
// Find infinite vertex and put it in k0
int k0 = 0;
int k1 = 1;
@ -298,13 +303,16 @@ interpolate_on_facet_vertices(const Bare_point& p, const Cell_handle& cell) cons
const FT& vb = cell->vertex(k2)->meshing_info();
const FT& vc = cell->vertex(k3)->meshing_info();
const Bare_point& a = wp2p(tr_.point(cell, k1));
const Bare_point& b = wp2p(tr_.point(cell, k2));
const Bare_point& c = wp2p(tr_.point(cell, k3));
const Weighted_point& wa = tr_.point(cell, k1);
const Weighted_point& wb = tr_.point(cell, k2);
const Weighted_point& wc = tr_.point(cell, k3);
const Bare_point& a = cp(wa);
const Bare_point& b = cp(wb);
const Bare_point& c = cp(wc);
const FT abp = area(a,b,p);
const FT acp = area(a,c,p);
const FT bcp = area(b,c,p);
const FT abp = area(a, b, p);
const FT acp = area(a, c, p);
const FT bcp = area(b, c, p);
CGAL_assertion(abp >= 0);
CGAL_assertion(acp >= 0);

View File

@ -49,6 +49,7 @@ class Odt_move
typedef typename Tr::Cell_handle Cell_handle;
typedef typename Tr::Bare_point Bare_point;
typedef typename Tr::Weighted_point Weighted_point;
typedef typename std::vector<Facet> Facet_vector;
typedef typename std::vector<Cell_handle> Cell_vector;
@ -73,10 +74,8 @@ public:
// Compute move
const Tr& tr = c3t3.triangulation();
typename Gt::Construct_vector_3 vector =
tr.geom_traits().construct_vector_3_object();
typename Gt::Construct_point_3 cp =
tr.geom_traits().construct_point_3_object();
typename Gt::Construct_point_3 cp = tr.geom_traits().construct_point_3_object();
typename Gt::Construct_vector_3 vector = tr.geom_traits().construct_vector_3_object();
Vector_3 move = CGAL::NULL_VECTOR;
FT sum_volume(0);
@ -95,8 +94,8 @@ public:
Bare_point circumcenter = tr.dual(cell);
// Compute move
const Bare_point& p = cp(tr.point(cell, cell->index(v)));
Vector_3 p_circum = vector(p, circumcenter);
const Weighted_point& p = tr.point(cell, cell->index(v));
Vector_3 p_circum = vector(cp(p), circumcenter);
FT volume = volume_quadrature(cell, tr, sizing_field);
move = move + p_circum * volume;
@ -126,13 +125,11 @@ private:
const Tr& tr,
const Sizing_field& sizing_field) const
{
typename Gt::Compute_volume_3 volume =
tr.geom_traits().compute_volume_3_object();
typename Gt::Construct_centroid_3 centroid =
tr.geom_traits().construct_centroid_3_object();
typename Gt::Construct_centroid_3 centroid = tr.geom_traits().construct_centroid_3_object();
typename Gt::Compute_volume_3 volume = tr.geom_traits().compute_volume_3_object();
Bare_point c = centroid(tr.tetrahedron(cell));
FT s = sizing_field(c,std::make_pair(cell,true));
FT s = sizing_field(c, std::make_pair(cell, true));
CGAL_assertion(!is_zero(s));
// Points of cell are positively oriented
@ -206,10 +203,8 @@ private:
// const Tr& tr,
// const Sizing_field& sizing_field) const
// {
// typename Gt::Compute_area_3 area =
// tr.geom_traits().compute_area_3_object();
// typename Gt::Construct_centroid_3 centroid =
// tr.geom_traits().construct_centroid_3_object();
// typename Gt::Compute_area_3 area = tr.geom_traits().compute_area_3_object();
// typename Gt::Construct_centroid_3 centroid = tr.geom_traits().construct_centroid_3_object();
//
// Bare_point c = centroid(tr.triangle(facet));
// FT s = sizing_field(c, facet.first->vertex(0));
@ -224,15 +219,14 @@ private:
// const Tr& tr,
// const Sizing_field& sizing_field) const
// {
// typename Gt::Construct_point_3 cp =
// tr.geom_traits().construct_point_3_object();
// typename Gt::Compute_squared_distance_3 sq_distance =
// tr.geom_traits().compute_squared_distance_3_object();
// typename Gt::Construct_midpoint_3 midpoint =
// tr.geom_traits().construct_midpoint_3_object();
// typename Gt::Construct_point_3 cp = tr.geom_traits().construct_point_3_object();
// typename Gt::Construct_midpoint_3 midpoint = tr.geom_traits().construct_midpoint_3_object();
// typename Gt::Compute_squared_distance_3 sq_distance = tr.geom_traits().compute_squared_distance_3_object();
//
// const Bare_point& p1 = cp(tr.point(cell, vertex_index_1));
// const Bare_point& p2 = cp(tr.point(cell, vertex_index_2));
// const Weighted_point& wp1 = tr.point(cell, vertex_index_1);
// const Weighted_point& wp2 = tr.point(cell, vertex_index_2);
// const Bare_point& p1 = cp(wp1);
// const Bare_point& p2 = cp(wp2);
//
// Bare_point c = midpoint(p1, p2);
// FT s = sizing_field(c, cell->vertex(vertex_index_1));
@ -265,10 +259,8 @@ private:
//
// Vector_3 normal_outside(const Facet& f, const C3T3& c3t3) const
// {
// typename Gt::Construct_point_3 cp =
// c3t3.triangulation().geom_traits().construct_point_3_object();
// typename Gt::Construct_normal_3 normal =
// c3t3.triangulation().geom_traits().construct_normal_3_object();
// typename Gt::Construct_point_3 cp = c3t3.triangulation().geom_traits().construct_point_3_object();
// typename Gt::Construct_normal_3 normal = c3t3.triangulation().geom_traits().construct_normal_3_object();
//
// const Cell_handle& cell = f.first;
// const int& i = f.second;
@ -285,11 +277,11 @@ private:
// if ( ! c3t3.is_in_complex(cell) )
// std::swap(k1,k2);
//
// const Bare_point& p1 = cp(c3t3.triangulation().point(cell, k1));
// const Bare_point& p2 = cp(c3t3.triangulation().point(cell, k2));
// const Bare_point& p3 = cp(c3t3.triangulation().point(cell, k3));
// const Weighted_point& wp1 = c3t3.triangulation().point(cell, k1);
// const Weighted_point& wp2 = c3t3.triangulation().point(cell, k2);
// const Weighted_point& wp3 = c3t3.triangulation().point(cell, k3);
//
// return normal(p1,p2,p3);
// return normal(cp(p1), cp(p2), cp(p3));
// }
};

View File

@ -109,11 +109,11 @@ struct Get_Is_facet_bad<Facet_criteria, true> {
#ifdef _DEBUG
int f1_current_erase_counter = CGAL::cpp11::get<0>(f).first->erase_counter();
int f1_saved_erase_counter = CGAL::cpp11::get<1>(f);
int f2_current_erase_counte7r = CGAL::cpp11::get<2>(f).first->erase_counter();
int f2_current_erase_counter = CGAL::cpp11::get<2>(f).first->erase_counter();
int f2_saved_erase_counter = CGAL::cpp11::get<3>(f);
//f1_current_erase_counter - f1_saved_erase_counter + f2_current_erase_counter - f2_saved_erase_counter == 1
/*if (f1_current_erase_counter - f1_saved_erase_counter + f2_current_erase_counter - f2_saved_erase_counter == 1)
if (f1_current_erase_counter - f1_saved_erase_counter + f2_current_erase_counter - f2_saved_erase_counter == 1)
{
#ifdef CGAL_LINKED_WITH_TBB
tbb::queuing_mutex mut;
@ -138,8 +138,8 @@ struct Get_Is_facet_bad<Facet_criteria, true> {
<< "}" << std::endl;
std::string s = sstr.str();
//std::cerr << s << std::endl;
}*/
std::cerr << s << std::endl;
}
#endif
return (CGAL::cpp11::get<0>(f).first->erase_counter() == CGAL::cpp11::get<1>(f)
&& CGAL::cpp11::get<2>(f).first->erase_counter() == CGAL::cpp11::get<3>(f) );
@ -297,12 +297,16 @@ public:
#ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS
const Cell_handle c = facet.first;
const int i = facet.second;
std::cerr << "Facet ("
std::cerr << "\n Facet ("
<< r_tr_.point(c, (i+1)&3) << " , "
<< r_tr_.point(c, (i+2)&3) << " , "
<< r_tr_.point(c, (i+3)&3) << ") : refinement point is "
<< get_facet_surface_center(facet) << std::endl;
#endif
if(CGAL::squared_distance(get_facet_surface_center(facet),
r_tr_.point(c, (i+3)&3)) < 1e-6)
exit(0);
CGAL_assertion (this->is_facet_on_surface(facet));
this->set_last_vertex_index(get_facet_surface_center_index(facet));
return get_facet_surface_center(facet);
@ -349,9 +353,9 @@ public:
{
std::stringstream sstr;
sstr << "Facet { " << std::endl
<< " " << *facet.first->vertex((facet.second+1)%4) << std::endl
<< " " << *facet.first->vertex((facet.second+2)%4) << std::endl
<< " " << *facet.first->vertex((facet.second+3)%4) << std::endl
<< " " << *facet.first->vertex((facet.second+1)%4) << std::endl
<< " " << *facet.first->vertex((facet.second+2)%4) << std::endl
<< " " << *facet.first->vertex((facet.second+3)%4) << std::endl
<< " 4th vertex in cell: " << *facet.first->vertex(facet.second) << std::endl
<< "}" << std::endl;
@ -973,9 +977,10 @@ scan_triangulation_impl()
{
// Cannot be const, see treat_new_facet signature
Facet facet = *facet_it;
/*std::cerr << "*" << *facet.first->vertex((facet.second+1)%4) << std::endl
<< " " << *facet.first->vertex((facet.second+2)%4) << std::endl
<< " " << *facet.first->vertex((facet.second+3)%4) << std::endl;*/
std::cerr << "TREAT FACET : " << std::endl
<< "*" << *facet.first->vertex((facet.second+1)%4) << std::endl
<< " " << *facet.first->vertex((facet.second+2)%4) << std::endl
<< " " << *facet.first->vertex((facet.second+3)%4) << std::endl;
this->treat_new_facet(facet);
}
}
@ -1388,10 +1393,10 @@ before_insertion_impl(const Facet& facet,
bool source_facet_is_in_conflict = false;
/*std::cerr << "before_insertion_impl:" << std::endl
<< "* " << *facet.first->vertex((facet.second+1)%4) << std::endl
<< " " << *facet.first->vertex((facet.second+2)%4) << std::endl
<< " " << *facet.first->vertex((facet.second+3)%4) << std::endl;*/
std::cerr << "before_insertion_impl:" << std::endl
<< "* " << *facet.first->vertex((facet.second+1)%4) << std::endl
<< " " << *facet.first->vertex((facet.second+2)%4) << std::endl
<< " " << *facet.first->vertex((facet.second+3)%4) << std::endl;
// Iterate on conflict zone facets
for (Facets_iterator facet_it = zone.internal_facets.begin();
@ -1544,11 +1549,12 @@ treat_new_facet(Facet& facet)
{
insert_bad_facet(facet, *is_facet_bad);
/*std::cerr << "INSERT BAD FACET : " << std::endl
<< "* " << *facet.first->vertex((facet.second+1)%4) << std::endl
<< " " << *facet.first->vertex((facet.second+2)%4) << std::endl
<< " " << *facet.first->vertex((facet.second+3)%4) << std::endl
<< " Quality=" << is_facet_bad->second << std::endl;*/
std::cerr << "INSERT BAD FACET : " << std::endl
<< "* " << *facet.first->vertex((facet.second+1)%4) << std::endl
<< " " << *facet.first->vertex((facet.second+2)%4) << std::endl
<< " " << *facet.first->vertex((facet.second+3)%4) << std::endl
<< " Surface center: " << get_facet_surface_center(facet) << std::endl
<< " Quality = " << is_facet_bad->second << std::endl;
}
}
else
@ -1741,9 +1747,9 @@ is_encroached_facet_refinable(Facet& facet) const
++wp_nb;
}
const Weighted_point& p1 = c->vertex(k1)->point();
const Weighted_point& p2 = c->vertex(k2)->point();
const Weighted_point& p3 = c->vertex(k3)->point();
const Weighted_point& p1 = r_tr_.point(c, k1);
const Weighted_point& p2 = r_tr_.point(c, k2);
const Weighted_point& p3 = r_tr_.point(c, k3);
const FT min_ratio (0.16); // (0.2*2)^2

View File

@ -61,29 +61,31 @@ class Refine_facets_manifold_base
Container_> Base;
public:
typedef Complex3InTriangulation3 C3t3;
typedef MeshDomain Mesh_domain;
typedef typename Tr::Weighted_point Weighted_point;
typedef typename Tr::Bare_point Bare_point;
typedef typename Tr::Facet Facet;
typedef typename Tr::Vertex_handle Vertex_handle;
typedef Complex3InTriangulation3 C3t3;
typedef MeshDomain Mesh_domain;
typedef typename Tr::Bare_point Bare_point;
typedef typename Tr::Weighted_point Weighted_point;
typedef typename Tr::Vertex_handle Vertex_handle;
typedef typename Tr::Facet Facet;
typedef typename Triangulation_mesher_level_traits_3<Tr>::Zone Zone;
protected:
typedef typename Tr::Geom_traits GT;
typedef typename GT::FT FT;
typedef typename GT::Construct_point_3 Construct_point_3;
typedef typename Tr::Edge Edge;
typedef typename Tr::Cell_handle Cell_handle;
typedef typename Tr::Geom_traits GT;
typedef typename GT::FT FT;
typedef typename Tr::Facet_circulator Tr_facet_circulator;
typedef std::pair<Vertex_handle, Vertex_handle> EdgeVV;
typedef typename Tr::Edge Edge;
typedef typename Tr::Cell_handle Cell_handle;
typedef typename Tr::Facet_circulator Tr_facet_circulator;
typedef std::pair<Vertex_handle, Vertex_handle> EdgeVV;
protected:
typedef ::boost::bimap< EdgeVV,
::boost::bimaps::multiset_of<int> > Bad_edges;
typedef typename Bad_edges::value_type Bad_edge;
::boost::bimaps::multiset_of<int> > Bad_edges;
typedef typename Bad_edges::value_type Bad_edge;
mutable Bad_edges m_bad_edges;
mutable std::set<Vertex_handle> m_bad_vertices;
@ -127,22 +129,21 @@ private:
}
FT compute_sq_distance_to_facet_center(const Facet& f,
const Vertex_handle v) const {
Construct_point_3 wp2p = this->r_tr_.geom_traits().construct_point_3_object();
const Bare_point& fcenter = f.first->get_facet_surface_center(f.second);
const Bare_point& vpoint = wp2p(v->point());
const Vertex_handle v) const
{
typename GT::Construct_point_3 cp = this->r_tr_.geom_traits().construct_point_3_object();
return
this->r_tr_.geom_traits().compute_squared_distance_3_object()(fcenter,
vpoint)
- v->point().weight();
const Bare_point& fcenter = f.first->get_facet_surface_center(f.second);
const Weighted_point& wp = this->r_tr_.point(v);
return this->r_tr_.min_squared_distance(fcenter, cp(wp)) - v->point().weight();
}
Facet
biggest_incident_facet_in_complex(const Vertex_handle v) const
{
#ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS
std::cerr << "Bad vertex: " << v->point() << std::endl;
std::cerr << "Bad vertex: " << this->r_tr_.point(v) << std::endl;
#endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS
std::vector<Facet> facets;
facets.reserve(64);
@ -201,8 +202,8 @@ private:
// use the list of incident facets in the complex
Vertex_handle fev = edge_to_edgevv(arete).first;
#ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS
std::cerr << "Bad edge: (" << fev->point()
<< ", " << arete.first->vertex(arete.third)->point()
std::cerr << "Bad edge: (" << this->r_tr_.point(fev)
<< ", " << this->r_tr_.point(arete.first->vertex(arete.third))
<< ")\n incident facets squared sizes:\n";
#endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS
Tr_facet_circulator fcirc = this->r_tr_.incident_facets(arete);
@ -290,7 +291,7 @@ private:
if(m_bad_vertices.erase(c->vertex(j)) > 0) {
#ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS
std::cerr << "m_bad_vertices.erase("
<< c->vertex(j)->point() << ")\n";
<< this->r_tr_.point(c, j) << ")\n";
#endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS
}
}
@ -345,18 +346,18 @@ public:
(this->r_c3t3_.face_status(*eit) == C3t3::BOUNDARY) ) )
{
#ifdef CGAL_LINKED_WITH_TBB
// Parallel
if (boost::is_convertible<Concurrency_tag, Parallel_tag>::value)
{
this->insert_bad_facet(biggest_incident_facet_in_complex(*eit),
typename Base::Quality());
} else
// Parallel
if (boost::is_convertible<Concurrency_tag, Parallel_tag>::value)
{
this->insert_bad_facet(biggest_incident_facet_in_complex(*eit),
typename Base::Quality());
} else
#endif // CGAL_LINKED_WITH_TBB
{ // Sequential
m_bad_edges.insert(Bad_edge(edge_to_edgevv(*eit),
(this->r_c3t3_.face_status(*eit) ==
C3t3::SINGULAR ? 0 : 1)));
}
{ // Sequential
m_bad_edges.insert(Bad_edge(edge_to_edgevv(*eit),
(this->r_c3t3_.face_status(*eit) ==
C3t3::SINGULAR ? 0 : 1)));
}
#ifdef CGAL_MESH_3_VERBOSE
++n;
#endif
@ -385,19 +386,19 @@ public:
if( this->r_c3t3_.face_status(vit) == C3t3::SINGULAR ) {
#ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS
std::cerr << "m_bad_vertices.insert("
<< vit->point() << ")\n";
<< this->r_tr_.point(vit) << ")\n";
#endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS
#ifdef CGAL_LINKED_WITH_TBB
// Parallel
if (boost::is_convertible<Concurrency_tag, Parallel_tag>::value)
{
this->insert_bad_facet(biggest_incident_facet_in_complex(vit),
typename Base::Quality());
} else
// Parallel
if (boost::is_convertible<Concurrency_tag, Parallel_tag>::value)
{
this->insert_bad_facet(biggest_incident_facet_in_complex(vit),
typename Base::Quality());
} else
#endif // CGAL_LINKED_WITH_TBB
{ // Sequential
m_bad_vertices.insert( vit );
}
{ // Sequential
m_bad_vertices.insert( vit );
}
#ifdef CGAL_MESH_3_VERBOSE
++n;
#endif
@ -410,7 +411,7 @@ public:
std::cerr << "Bad vertices queue:\n";
BOOST_FOREACH(Vertex_handle v2, m_bad_vertices)
{
std::cerr << v2->point() << std::endl;
std::cerr << this->r_tr_.point(v2) << std::endl;
}
CGAL::dump_c3t3(this->r_c3t3_, "dump-at-scan-vertices");
# endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS
@ -456,9 +457,9 @@ public:
#ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS
const EdgeVV& edgevv = m_bad_edges.right.begin()->second;
std::cerr << "Bad edge "
<< edgevv.first->point()
<< this->r_tr_.point(edgevv.first)
<< " - "
<< edgevv.second->point()
<< this->r_tr_.point(edgevv.second)
<< "\n";
#endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS
return biggest_incident_facet_in_complex(first_bad_edge);
@ -469,9 +470,9 @@ public:
std::cerr << "Bad vertices queue:\n";
BOOST_FOREACH(Vertex_handle v2, m_bad_vertices)
{
std::cerr << v2->point() << std::endl;
std::cerr << this->r_tr_.point(v2) << std::endl;
}
std::cerr << "Bad vertex " << v->point() << "\n";
std::cerr << "Bad vertex " << this->r_tr_.point(v) << "\n";
#endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS
CGAL_assertion(this->r_c3t3_.has_incident_facets_in_complex(v));
if(this->r_c3t3_.face_status(v) != C3t3::SINGULAR) {
@ -541,27 +542,27 @@ public:
)
{
#ifdef CGAL_LINKED_WITH_TBB
// Parallel
if (boost::is_convertible<Concurrency_tag, Parallel_tag>::value)
{
this->insert_bad_facet(biggest_incident_facet_in_complex(edge),
typename Base::Quality());
} else
// Parallel
if (boost::is_convertible<Concurrency_tag, Parallel_tag>::value)
{
this->insert_bad_facet(biggest_incident_facet_in_complex(edge),
typename Base::Quality());
} else
#endif // CGAL_LINKED_WITH_TBB
{ // Sequential
m_bad_edges.insert(Bad_edge(edge_to_edgevv(edge),
(this->r_c3t3_.face_status(edge) ==
C3t3::SINGULAR ? 0 : 1)));
}
{ // Sequential
m_bad_edges.insert(Bad_edge(edge_to_edgevv(edge),
(this->r_c3t3_.face_status(edge) ==
C3t3::SINGULAR ? 0 : 1)));
}
}
else {
#ifdef CGAL_LINKED_WITH_TBB
// Sequential only
if (!boost::is_convertible<Concurrency_tag, Parallel_tag>::value)
// Sequential only
if (!boost::is_convertible<Concurrency_tag, Parallel_tag>::value)
#endif // CGAL_LINKED_WITH_TBB
{
m_bad_edges.left.erase( edge_to_edgevv(edge) ); // @TODO: pourquoi?!
}
{
m_bad_edges.left.erase( edge_to_edgevv(edge) ); // @TODO: pourquoi?!
}
}
}
}
@ -590,19 +591,19 @@ public:
{
#ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS
std::cerr << "m_bad_vertices.insert("
<< (*vit)->point() << ")\n";
<< this->r_tr_.point(*vit) << ")\n";
#endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS
#ifdef CGAL_LINKED_WITH_TBB
// Parallel
if (boost::is_convertible<Concurrency_tag, Parallel_tag>::value)
{
this->insert_bad_facet(biggest_incident_facet_in_complex(*vit),
typename Base::Quality());
} else
// Parallel
if (boost::is_convertible<Concurrency_tag, Parallel_tag>::value)
{
this->insert_bad_facet(biggest_incident_facet_in_complex(*vit),
typename Base::Quality());
} else
#endif // CGAL_LINKED_WITH_TBB
{ // Sequential
m_bad_vertices.insert(*vit);
}
{ // Sequential
m_bad_vertices.insert(*vit);
}
}
}
@ -613,19 +614,19 @@ public:
{
#ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS
std::cerr << "m_bad_vertices.insert("
<< v->point() << ")\n";
<< this->r_tr_.point(v) << ")\n";
#endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS
#ifdef CGAL_LINKED_WITH_TBB
// Parallel
if (boost::is_convertible<Concurrency_tag, Parallel_tag>::value)
{
this->insert_bad_facet(biggest_incident_facet_in_complex(v),
typename Base::Quality());
} else
// Parallel
if (boost::is_convertible<Concurrency_tag, Parallel_tag>::value)
{
this->insert_bad_facet(biggest_incident_facet_in_complex(v),
typename Base::Quality());
} else
#endif // CGAL_LINKED_WITH_TBB
{ // Sequential
m_bad_vertices.insert(v);
}
{ // Sequential
m_bad_vertices.insert(v);
}
}
}

View File

@ -472,13 +472,15 @@ class Sliver_perturber
typedef Sliver_perturber_base<
typename C3T3::Triangulation, Concurrency_tag> Base;
typedef typename C3T3::Triangulation Tr;
typedef typename Tr::Geom_traits Gt;
typedef typename C3T3::Triangulation Tr;
typedef typename Tr::Geom_traits Gt;
typedef typename Tr::Cell_handle Cell_handle;
typedef typename Base::Vertex_handle Vertex_handle;
typedef typename Tr::Vertex Vertex;
typedef typename MeshDomain::Point_3 Point_3;
typedef typename Tr::Bare_point Bare_point;
typedef typename Tr::Weighted_point Weighted_point;
typedef typename std::vector<Cell_handle> Cell_vector;
typedef typename std::vector<Vertex_handle> Vertex_vector;
@ -507,7 +509,7 @@ private:
typedef PVertex_<FT,
Vertex_handle,
Point_3,
Bare_point,
SliverCriterion,
Perturbation,
Concurrency_tag> PVertex;
@ -1298,7 +1300,7 @@ perturb_vertex( PVertex pv
, bool *could_lock_zone
) const
{
typename Gt::Construct_point_3 wp2p = tr_.geom_traits().construct_point_3_object();
typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object();
#ifdef CGAL_CONCURRENT_MESH_3_PROFILING
static Profile_branch_counter_3 bcounter(
@ -1313,9 +1315,10 @@ perturb_vertex( PVertex pv
return;
}
Point_3 p = wp2p(pv.vertex()->point());
// <periodic> and what's happening below... ?
Point_3 p = cp(pv.vertex()->point());
if (!helper_.try_lock_point_no_spin(pv.vertex()->point()) ||
! tr_.geom_traits().equal_3_object()(p, wp2p(pv.vertex()->point())))
! tr_.geom_traits().equal_3_object()(p, cp(pv.vertex()->point())))
{
#ifdef CGAL_CONCURRENT_MESH_3_PROFILING
bcounter.increment_branch_2(); // THIS is an early withdrawal!

View File

@ -99,34 +99,34 @@ namespace Mesh_3 {
// Then, its operator() takes an other vertex handle v2 as input, and
// returns the distance d(v1, v2).
// It is used in Slivers_exuder, to constructor a transform iterator.
template <typename Gt, typename Vertex_handle>
template <typename Tr, typename Vertex_handle>
class Min_distance_from_v :
public CGAL::unary_function<Vertex_handle, void>
public CGAL::unary_function<Vertex_handle, void>
{
const Vertex_handle * v;
const Gt& gt;
const Tr& tr;
double & dist;
public:
Min_distance_from_v(const Vertex_handle& vh,
double& dist,
const Gt& geom_traits = Gt())
: v(&vh), gt(geom_traits), dist(dist)
const Tr& tr)
: v(&vh), tr(tr), dist(dist)
{
}
void
operator()(const Vertex_handle& vh) const
{
typedef typename Gt::Compute_squared_distance_3 Compute_squared_distance_3;
Compute_squared_distance_3 distance = gt.compute_squared_distance_3_object();
typedef typename Tr::Geom_traits::Construct_point_3 Construct_point_3;
typedef typename Tr::Weighted_point Weighted_point;
typedef typename Gt::Construct_point_3 Construct_point_3;
Construct_point_3 wp2p = gt.construct_point_3_object();
Construct_point_3 cp = tr.geom_traits().construct_point_3_object();
// <periodic...>
const double d = CGAL::to_double(distance(wp2p((*v)->point()),
wp2p(vh->point())));
const Weighted_point& wpv = tr.point(*v);
const Weighted_point& wpvh = tr.point(vh);
const double d = CGAL::to_double(tr.min_squared_distance(cp(wpv), cp(wpvh)));
if(d < dist){
dist = d;
}
@ -645,7 +645,7 @@ private:
{
double dist = (std::numeric_limits<double>::max)();
details::Min_distance_from_v<Gt, Vertex_handle> min_distance_from_v(vh, dist, tr_.geom_traits());
details::Min_distance_from_v<Tr, Vertex_handle> min_distance_from_v(vh, dist, tr_);
tr_.adjacent_vertices(vh, boost::make_function_output_iterator(min_distance_from_v));
@ -1126,8 +1126,10 @@ pump_vertex(const Vertex_handle& pumped_vertex,
// If best_weight < pumped_vertex weight, nothing to do
if ( best_weight > pumped_vertex->point().weight() )
{
typename Gt::Construct_point_3 wp2p = tr_.geom_traits().construct_point_3_object();
Weighted_point wp(wp2p(pumped_vertex->point()), best_weight);
typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object();
const Weighted_point& pwp = tr_.point(pumped_vertex);
Weighted_point wp(cp(pwp), best_weight);
// Insert weighted point into mesh
// note it can fail if the mesh is non-manifold at pumped_vertex
@ -1205,6 +1207,8 @@ expand_prestar(const Cell_handle& cell_to_add,
Pre_star& pre_star,
Sliver_values& criterion_values) const
{
typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object();
// Delete first facet of pre_star
Facet start_facet = pre_star.front()->second;
CGAL_assertion(tr_.mirror_facet(start_facet).first == cell_to_add);
@ -1278,14 +1282,15 @@ expand_prestar(const Cell_handle& cell_to_add,
// Update ratio (ratio is needed for cells of complex only)
if ( c3t3_.is_in_complex(cell_to_add) )
{
typename Gt::Construct_point_3 wp2p = tr_.geom_traits().construct_point_3_object();
Tetrahedron_3 tet(wp2p(pumped_vertex->point()),
wp2p(cell_to_add->vertex((i+1)&3)->point()),
wp2p(cell_to_add->vertex((i+2)&3)->point()),
wp2p(cell_to_add->vertex((i+3)&3)->point()));
const Weighted_point& pwp = tr_.point(pumped_vertex);
const Weighted_point& fwp1 = tr_.point(cell_to_add, (i+1)&3);
const Weighted_point& fwp2 = tr_.point(cell_to_add, (i+2)&3);
const Weighted_point& fwp3 = tr_.point(cell_to_add, (i+3)&3);
const Tetrahedron_3 tet(cp(pwp), cp(fwp1), cp(fwp2), cp(fwp3));
double new_value = sliver_criteria_(tet);
criterion_values.insert(std::make_pair(current_facet,new_value));
criterion_values.insert(std::make_pair(current_facet, new_value));
}
}
}
@ -1369,7 +1374,7 @@ get_best_weight(const Vertex_handle& v, bool *could_lock_zone) const
#ifdef CGAL_MESH_3_DEBUG_SLIVERS_EXUDER
if ( best_weight > v->point().weight() )
{
Weighted_point wp(v->point(), best_weight);
Weighted_point wp(tr_.point(v), best_weight);
check_pre_star(pre_star_copy, wp, v);
check_ratios(ratios_copy, wp, v);
}
@ -1794,7 +1799,7 @@ check_ratios(const Sliver_values& criterion_values,
Facet_vector internal_facets;
Facet_vector boundary_facets;
typename Gt::Construct_point_3 wp2p = tr_.geom_traits().construct_point_3_object();
typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object();
tr_.find_conflicts(wp,
vh->cell(),
@ -1821,10 +1826,13 @@ check_ratios(const Sliver_values& criterion_values,
continue;
int k = it->second;
Tetrahedron_3 tet(wp2p(vh->point()),
wp2p(it->first->vertex((k+1)&3)->point()),
wp2p(it->first->vertex((k+2)&3)->point()),
wp2p(it->first->vertex((k+3)&3)->point()));
const Weighted_point& vhwp = tr_.point(vh);
const Weighted_point& fwp1 = tr_.point(it->first, (k+1)&3));
const Weighted_point& fwp2 = tr_.point(it->first, (k+2)&3));
const Weighted_point& fwp3 = tr_.point(it->first, (k+3)&3));
Tetrahedron_3 tet(cp(vhwp), cp(fwp1), cp(fwp2), cp(fwp3));
double ratio = sliver_criteria_(tet);
expected_ratios.push_back(ratio);

View File

@ -172,7 +172,7 @@ no_topological_change(Tr& tr,
tr.geom_traits().construct_opposite_vector_3_object();
bool np = true;
const Weighted_point fp = v0->point();
const Weighted_point fp = tr.point(v0);
//#define CGAL_PERIODIC_BACKUP_CHECK
//#define CGAL_PERIODIC_SIDE_OF_DEBUG
@ -528,12 +528,16 @@ inside_protecting_balls(const Tr& tr,
{
Vertex_handle nv = tr.nearest_power_vertex(p, v->cell());
if(nv->point().weight() > 0)
return tr.geom_traits().compare_squared_distance_3_object()(
p, nv->point(), nv->point().weight()) != CGAL::LARGER;
{
typename Tr::Geom_traits::Construct_point_3 cp = tr.geom_traits().construct_point_3_object();
const Weighted_point& nvwp = tr.point(nv);
return (tr.min_squared_distance(p, cp(nvwp)) <= nv->point().weight());
}
return false;
}
/// This function well_oriented is called by no_topological_change after the
/// position of the vertex has been (tentatively) modified.
template<typename Tr>
@ -544,7 +548,7 @@ well_oriented(const Tr& tr,
{
typedef typename Tr::Geom_traits Gt;
typename Gt::Orientation_3 orientation = tr.geom_traits().orientation_3_object();
typename Gt::Construct_point_3 wp2p = tr.geom_traits().construct_point_3_object();
typename Gt::Construct_point_3 cp = tr.geom_traits().construct_point_3_object();
typename Cell_vector::const_iterator it = cells_tos.begin();
for( ; it != cells_tos.end() ; ++it)
@ -555,17 +559,25 @@ well_oriented(const Tr& tr,
int iv = c->index(tr.infinite_vertex());
Cell_handle cj = c->neighbor(iv);
int mj = tr.mirror_index(c, iv);
if(orientation(wp2p(tr.point(cj, mj)),
wp2p(tr.point(c, (iv+1)&3)),
wp2p(tr.point(c, (iv+2)&3)),
wp2p(tr.point(c, (iv+3)&3))) != CGAL::NEGATIVE)
const Weighted_point& mjwp = tr.point(cj, mj);
const Weighted_point& fwp1 = tr.point(c, (iv+1)&3);
const Weighted_point& fwp2 = tr.point(c, (iv+2)&3);
const Weighted_point& fwp3 = tr.point(c, (iv+3)&3);
if(orientation(cp(mjwp), cp(fwp1), cp(fwp2), cp(fwp3)) != CGAL::NEGATIVE)
return false;
}
else if(orientation(wp2p(tr.point(c, 0)),
wp2p(tr.point(c, 1)),
wp2p(tr.point(c, 2)),
wp2p(tr.point(c, 3))) != CGAL::POSITIVE)
else
{
const Weighted_point& cwp0 = tr.point(c, 0);
const Weighted_point& cwp1 = tr.point(c, 1);
const Weighted_point& cwp2 = tr.point(c, 2);
const Weighted_point& cwp3 = tr.point(c, 3);
if(orientation(cp(cwp0), cp(cwp1), cp(cwp2), cp(cwp3)) != CGAL::POSITIVE)
return false;
}
}
return true;
}
@ -582,7 +594,7 @@ well_oriented(const Tr& tr,
{
typedef typename Tr::Geom_traits Gt;
typename Gt::Orientation_3 orientation = tr.geom_traits().orientation_3_object();
typename Gt::Construct_point_3 wp2p = tr.geom_traits().construct_point_3_object();
typename Gt::Construct_point_3 cp = tr.geom_traits().construct_point_3_object();
typename Cell_vector::const_iterator it = cells_tos.begin();
for( ; it != cells_tos.end() ; ++it)
@ -593,16 +605,16 @@ well_oriented(const Tr& tr,
int iv = c->index(tr.infinite_vertex());
Cell_handle cj = c->neighbor(iv);
int mj = tr.mirror_index(c, iv);
if(orientation(wp2p(pg(cj->vertex(mj))),
wp2p(pg(c->vertex((iv+1)&3))),
wp2p(pg(c->vertex((iv+2)&3))),
wp2p(pg(c->vertex((iv+3)&3)))) != CGAL::NEGATIVE)
if(orientation(cp(pg(cj->vertex(mj))),
cp(pg(c->vertex((iv+1)&3))),
cp(pg(c->vertex((iv+2)&3))),
cp(pg(c->vertex((iv+3)&3)))) != CGAL::NEGATIVE)
return false;
}
else if(orientation(wp2p(pg(c->vertex(0))),
wp2p(pg(c->vertex(1))),
wp2p(pg(c->vertex(2))),
wp2p(pg(c->vertex(3)))) != CGAL::POSITIVE)
else if(orientation(cp(pg(c->vertex(0))),
cp(pg(c->vertex(1))),
cp(pg(c->vertex(2))),
cp(pg(c->vertex(3)))) != CGAL::POSITIVE)
return false;
}
return true;

View File

@ -40,7 +40,6 @@ namespace CGAL {
namespace Mesh_3
{
/**
* @class Triangulation_sizing_field
@ -49,20 +48,21 @@ template <typename Tr>
class Triangulation_sizing_field
{
// Types
typedef typename Tr::Geom_traits Gt;
typedef typename Tr::Weighted_point Weighted_point;
typedef typename Gt::FT FT;
typedef typename Tr::Geom_traits Gt;
typedef typename Tr::Bare_point Bare_point;
typedef typename Tr::Weighted_point Weighted_point;
typedef typename Gt::FT FT;
typedef Triangulation_vertex_base_with_info_3<FT, Gt> Vb;
typedef Triangulation_cell_base_3<Gt> Cb;
typedef Triangulation_data_structure_3<Vb, Cb> Tds;
typedef Regular_triangulation_3<Gt,Tds> Compact_triangulation;
typedef Compact_triangulation Ctr;
typedef typename Tr::Vertex_handle Vertex_handle;
typedef typename Tr::Cell_handle Cell_handle;
typedef typename Ctr::Cell_handle CCell_handle;
typedef typename Ctr::Vertex_handle CVertex_handle;
typedef typename Tr::Vertex_handle Vertex_handle;
typedef typename Tr::Cell_handle Cell_handle;
typedef typename Ctr::Cell_handle CCell_handle;
typedef typename Ctr::Vertex_handle CVertex_handle;
public:
// Vertices of mesh triangulation do not need to be updated
@ -116,8 +116,9 @@ private:
* Used by boost transform iterator
*/
struct Extract_point :
public CGAL::unary_function<typename Tr::Vertex,Weighted_point>
public CGAL::unary_function<typename Tr::Vertex, Weighted_point>
{
// <periodic>
Weighted_point operator()(const typename Tr::Vertex& v) const { return v.point(); }
};
@ -152,7 +153,7 @@ fill(const std::map<Weighted_point, FT>& value_map)
++ vit )
{
typename std::map<Weighted_point, FT>::const_iterator find_result =
value_map.find(vit->point());
value_map.find(ctr_.point(vit));
if ( find_result != value_map.end() )
vit->info() = find_result->second;
@ -183,20 +184,24 @@ typename Triangulation_sizing_field<Tr>::FT
Triangulation_sizing_field<Tr>::
interpolate_on_cell_vertices(const Weighted_point& p, const CCell_handle& cell) const
{
typename Gt::Compute_volume_3 volume =
ctr_.geom_traits().compute_volume_3_object();
typename Gt::Construct_point_3 cp = ctr_.geom_traits().construct_point_3_object();
typename Gt::Compute_volume_3 volume = ctr_.geom_traits().compute_volume_3_object();
// Interpolate value using tet vertices values
const FT& va = cell->vertex(0)->info();
const FT& vb = cell->vertex(1)->info();
const FT& vc = cell->vertex(2)->info();
const FT& vd = cell->vertex(3)->info();
const Weighted_point& a = cell->vertex(0)->point();
const Weighted_point& b = cell->vertex(1)->point();
const Weighted_point& c = cell->vertex(2)->point();
const Weighted_point& d = cell->vertex(3)->point();
const Weighted_point& wa = ctr_.point(cell, 0);
const Weighted_point& wb = ctr_.point(cell, 1);
const Weighted_point& wc = ctr_.point(cell, 2);
const Weighted_point& wd = ctr_.point(cell, 3);
const Bare_point& a = cp(wa);
const Bare_point& b = cp(wb);
const Bare_point& c = cp(wc);
const Bare_point& d = cp(wd);
const FT abcp = CGAL::abs(volume(a,b,c,p));
const FT abdp = CGAL::abs(volume(a,b,d,p));
const FT acdp = CGAL::abs(volume(a,c,d,p));
@ -215,8 +220,8 @@ typename Triangulation_sizing_field<Tr>::FT
Triangulation_sizing_field<Tr>::
interpolate_on_facet_vertices(const Weighted_point& p, const CCell_handle& cell) const
{
typename Gt::Compute_area_3 area =
ctr_.geom_traits().compute_area_3_object();
typename Gt::Construct_point_3 cp = ctr_.geom_traits().construct_point_3_object();
typename Gt::Compute_area_3 area = ctr_.geom_traits().compute_area_3_object();
// Find infinite vertex and put it in k0
int k0 = 0;
@ -236,10 +241,13 @@ interpolate_on_facet_vertices(const Weighted_point& p, const CCell_handle& cell)
const FT& vb = cell->vertex(k2)->info();
const FT& vc = cell->vertex(k3)->info();
const Weighted_point& a = cell->vertex(k1)->point();
const Weighted_point& b = cell->vertex(k2)->point();
const Weighted_point& c = cell->vertex(k3)->point();
const Weighted_point& wa = ctr_.point(cell, k1);
const Weighted_point& wb = ctr_.point(cell, k2);
const Weighted_point& wc = ctr_.point(cell, k3);
const Bare_point& a = cp(wa);
const Bare_point& b = cp(wb);
const Bare_point& c = cp(wc);
const FT abp = area(a,b,p);
const FT acp = area(a,c,p);
const FT bcp = area(b,c,p);

View File

@ -84,21 +84,21 @@ void initialize_triangulation_from_labeled_image(C3T3& c3t3,
)
{
typedef typename C3T3::Triangulation Tr;
typedef typename Tr::Geom_traits Gt;
typedef typename Tr::Weighted_point Weighted_point;
typedef typename Tr::Bare_point Bare_point;
typedef typename Tr::Segment Segment_3;
typedef typename Tr::Geom_traits::Vector_3 Vector_3;
typedef typename Tr::Vertex_handle Vertex_handle;
typedef typename Tr::Cell_handle Cell_handle;
typedef typename Gt::Vector_3 Vector_3;
typedef MeshDomain Mesh_domain;
Tr& tr = c3t3.triangulation();
typename Tr::Geom_traits::Construct_point_3 wp2p =
tr.geom_traits().construct_point_3_object();
typename Tr::Geom_traits::Construct_weighted_point_3 p2wp =
tr.geom_traits().construct_weighted_point_3_object();
typename Gt::Construct_point_3 cp = tr.geom_traits().construct_point_3_object();
typename Gt::Construct_weighted_point_3 cwp = tr.geom_traits().construct_weighted_point_3_object();
if(protect_features) {
init_tr_from_labeled_image_call_init_features
@ -155,7 +155,8 @@ void initialize_triangulation_from_labeled_image(C3T3& c3t3,
construct_intersection(Segment_3(it->first, test));
if (CGAL::cpp11::get<2>(intersect) != 0)
{
Weighted_point pi = p2wp(CGAL::cpp11::get<0>(intersect));
const Bare_point& bpi = CGAL::cpp11::get<0>(intersect);
Weighted_point pi = cwp(bpi);
// This would cause trouble to optimizers
// check pi will not be hidden
@ -201,8 +202,9 @@ void initialize_triangulation_from_labeled_image(C3T3& c3t3,
{
if (cv->point().weight() == 0.)
continue;
if (CGAL::compare_squared_distance(pi.point(), wp2p(cv->point()), cv->point().weight())
!= CGAL::LARGER)
const Weighted_point& cvwp = tr.point(cv);
if (tr.min_squared_distance(bpi, cp(cvwp) <= cv->point().weight()))
{
pi_inside_protecting_sphere = true;
break;

View File

@ -79,20 +79,26 @@ protected:
virtual Badness do_is_bad(const Tr& tr, const Cell_handle& ch) const
{
typedef typename Tr::Bare_point Bare_point;
typedef typename Tr::Geom_traits Geom_traits;
typedef typename Tr::Bare_point Bare_point;
typedef typename Tr::Weighted_point Weighted_point;
typedef typename Geom_traits::Compute_squared_distance_3 Distance;
typedef typename Geom_traits::Compute_squared_radius_3 Radius;
typedef typename Geom_traits::Construct_point_3 Construct_point_3;
Distance distance = tr.geom_traits().compute_squared_distance_3_object();
Radius sq_radius = tr.geom_traits().compute_squared_radius_3_object();
Construct_point_3 wp2p = tr.geom_traits().construct_point_3_object();
Construct_point_3 cp = tr.geom_traits().construct_point_3_object();
const Bare_point& p = wp2p(tr.point(ch, 0));
const Bare_point& q = wp2p(tr.point(ch, 1));
const Bare_point& r = wp2p(tr.point(ch, 2));
const Bare_point& s = wp2p(tr.point(ch, 3));
const Weighted_point& wp = tr.point(ch, 0);
const Weighted_point& wq = tr.point(ch, 1);
const Weighted_point& wr = tr.point(ch, 2);
const Weighted_point& ws = tr.point(ch, 3);
const Bare_point& p = cp(wp);
const Bare_point& q = cp(wq);
const Bare_point& r = cp(wr);
const Bare_point& s = cp(ws);
const FT size = sq_radius(p, q, r, s);
@ -164,18 +170,23 @@ protected:
virtual Badness do_is_bad(const Tr& tr, const Cell_handle& ch) const
{
typedef typename Tr::Bare_point Bare_point;
typedef typename Tr::Geom_traits Geom_traits;
typedef typename Tr::Geom_traits Geom_traits;
typedef typename Tr::Bare_point Bare_point;
typedef typename Tr::Weighted_point Weighted_point;
typedef typename Geom_traits::Compute_squared_radius_3 Radius;
typedef typename Geom_traits::Construct_point_3 Construct_point_3;
Radius sq_radius = tr.geom_traits().compute_squared_radius_3_object();
Construct_point_3 wp2p = tr.geom_traits().construct_point_3_object();
Construct_point_3 cp = tr.geom_traits().construct_point_3_object();
const Bare_point& p = wp2p(tr.point(ch, 0));
const Bare_point& q = wp2p(tr.point(ch, 1));
const Bare_point& r = wp2p(tr.point(ch, 2));
const Bare_point& s = wp2p(tr.point(ch, 3));
const Weighted_point& wp = tr.point(ch, 0);
const Weighted_point& wq = tr.point(ch, 1);
const Weighted_point& wr = tr.point(ch, 2);
const Weighted_point& ws = tr.point(ch, 3);
const Bare_point& p = cp(wp);
const Bare_point& q = cp(wq);
const Bare_point& r = cp(wr);
const Bare_point& s = cp(ws);
const FT size = sq_radius(p, q, r, s);
@ -234,18 +245,23 @@ protected:
virtual Badness do_is_bad(const Tr& tr, const Cell_handle& ch) const
{
typedef typename Tr::Bare_point Bare_point;
typedef typename Tr::Geom_traits Geom_traits;
typedef typename Tr::Geom_traits Geom_traits;
typedef typename Tr::Bare_point Bare_point;
typedef typename Tr::Weighted_point Weighted_point;
typedef typename Geom_traits::Compute_squared_radius_3 Radius;
typedef typename Geom_traits::Construct_point_3 Construct_point_3;
Radius sq_radius = tr.geom_traits().compute_squared_radius_3_object();
Construct_point_3 wp2p = tr.geom_traits().construct_point_3_object();
Construct_point_3 cp = tr.geom_traits().construct_point_3_object();
const Bare_point& p = wp2p(tr.point(ch, 0));
const Bare_point& q = wp2p(tr.point(ch, 1));
const Bare_point& r = wp2p(tr.point(ch, 2));
const Bare_point& s = wp2p(tr.point(ch, 3));
const Weighted_point& wp = tr.point(ch, 0);
const Weighted_point& wq = tr.point(ch, 1);
const Weighted_point& wr = tr.point(ch, 2);
const Weighted_point& ws = tr.point(ch, 3);
const Bare_point& p = cp(wp);
const Bare_point& q = cp(wq);
const Bare_point& r = cp(wr);
const Bare_point& s = cp(ws);
const FT size = sq_radius(p, q, r, s);
const FT sq_bound = CGAL::square( size_(tr.dual(ch), 3,

View File

@ -110,8 +110,9 @@ protected:
CGAL_assertion (f.first->is_facet_on_surface(f.second));
CGAL_assertion (B_ != 0);
typedef typename Tr::Geom_traits Gt;
typedef typename Tr::Bare_point Bare_point;
typedef typename Tr::Geom_traits Gt;
typedef typename Tr::Bare_point Bare_point;
typedef typename Tr::Weighted_point Weighted_point;
typedef typename Gt::Compute_squared_area_3 Area;
typedef typename Gt::Compute_squared_distance_3 Distance;
@ -120,12 +121,15 @@ protected:
Area area = tr.geom_traits().compute_squared_area_3_object();
Distance distance = tr.geom_traits().compute_squared_distance_3_object();
Construct_point_3 wp2p = tr.geom_traits().construct_point_3_object();
Construct_point_3 cp = tr.geom_traits().construct_point_3_object();
Construct_triangle_3 triangle = tr.geom_traits().construct_triangle_3_object();
const Bare_point& p1 = wp2p(tr.point(f.first, (f.second+1)&3));
const Bare_point& p2 = wp2p(tr.point(f.first, (f.second+2)&3));
const Bare_point& p3 = wp2p(tr.point(f.first, (f.second+3)&3));
const Weighted_point& wp1 = tr.point(f.first, (f.second+1)&3);
const Weighted_point& wp2 = tr.point(f.first, (f.second+2)&3);
const Weighted_point& wp3 = tr.point(f.first, (f.second+3)&3);
const Bare_point& p1 = cp(wp1);
const Bare_point& p2 = cp(wp2);
const Bare_point& p3 = cp(wp3);
const FT triangle_area = area(triangle(p1,p2,p3));
const FT d12 = distance(p1,p2);
@ -351,12 +355,13 @@ protected:
CGAL_assertion (f.first->is_facet_on_surface(f.second));
typedef typename Tr::Geom_traits Gt;
typedef typename Tr::Bare_point Bare_point;
typedef typename Tr::Bare_point Bare_point;
typedef typename Tr::Weighted_point Weighted_point;
typename Gt::Construct_point_3 wp2p =
tr.geom_traits().construct_point_3_object();
typename Gt::Construct_point_3 cp = tr.geom_traits().construct_point_3_object();
const Bare_point& p1 = wp2p(tr.point(f.first, (f.second+1)&3));
const Weighted_point& wp1 = tr.point(f.first, (f.second+1)&3);
const Bare_point& p1 = cp(wp1);
const Bare_point& ball_center = f.first->get_facet_surface_center(f.second);
const Index& index = f.first->get_facet_surface_center_index(f.second);
@ -420,13 +425,14 @@ protected:
CGAL_assertion (f.first->is_facet_on_surface(f.second));
CGAL_assertion (B_ != 0);
typedef typename Tr::Geom_traits Gt;
typedef typename Tr::Bare_point Bare_point;
typedef typename Tr::Geom_traits Gt;
typedef typename Tr::Bare_point Bare_point;
typedef typename Tr::Weighted_point Weighted_point;
typename Gt::Construct_point_3 wp2p =
tr.geom_traits().construct_point_3_object();
typename Gt::Construct_point_3 cp = tr.geom_traits().construct_point_3_object();
const Bare_point& p1 = wp2p(tr.point(f.first, (f.second+1)&3));
const Weighted_point& wp1 = tr.point(f.first, (f.second+1)&3);
const Bare_point p1 = cp(wp1);
const Bare_point& ball_center = f.first->get_facet_surface_center(f.second);
const FT sq_radius = tr.min_squared_distance(p1, ball_center);

View File

@ -123,16 +123,19 @@ namespace Mesh_3 {
edge_sq_length(const typename Tr::Edge& e,
const Tr& tr)
{
typedef typename Tr::Geom_traits Gt;
typedef typename Tr::Bare_point Bare_point;
typedef typename Tr::Geom_traits Gt;
typedef typename Tr::Bare_point Bare_point;
typedef typename Tr::Weighted_point Weighted_point;
typename Gt::Compute_squared_distance_3 sq_distance =
tr.geom_traits().compute_squared_distance_3_object();
typename Gt::Construct_point_3 wp2p =
typename Gt::Construct_point_3 cp =
tr.geom_traits().construct_point_3_object();
const Bare_point& p = wp2p(tr.point(e.first, e.second));
const Bare_point& q = wp2p(tr.point(e.first, e.third));
const Weighted_point& wp = tr.point(e.first, e.second);
const Weighted_point& wq = tr.point(e.first, e.third);
const Bare_point& p = cp(wp);
const Bare_point& q = cp(wq);
return sq_distance(p,q);
}
@ -182,10 +185,11 @@ template <typename C3T3, typename MeshDomain, typename SliverCriterion>
class Abstract_perturbation
{
protected:
typedef typename C3T3::Vertex_handle Vertex_handle;
typedef typename C3T3::Cell_handle Cell_handle;
typedef typename C3T3::Triangulation::Geom_traits::FT FT;
typedef typename C3T3::Vertex_handle Vertex_handle;
typedef typename C3T3::Cell_handle Cell_handle;
typedef typename C3T3::Triangulation::Geom_traits::FT FT;
public:
/**
* @brief constructor
@ -393,12 +397,19 @@ class Gradient_based_perturbation
: public Abstract_perturbation<C3T3,MeshDomain,SliverCriterion>
{
protected:
typedef Abstract_perturbation<C3T3,MeshDomain,SliverCriterion> Base;
typedef typename Base::Vertex_handle Vertex_handle;
typedef typename Base::Cell_handle Cell_handle;
typedef typename Base::FT FT;
typedef typename C3T3::Triangulation::Geom_traits::Vector_3 Vector_3;
typedef Abstract_perturbation<C3T3, MeshDomain, SliverCriterion> Base;
typedef typename Base::Vertex_handle Vertex_handle;
typedef typename Base::Cell_handle Cell_handle;
typedef typename C3T3::Triangulation Tr;
typedef typename Tr::Geom_traits Gt;
typedef typename Gt::FT FT;
typedef typename Gt::Vector_3 Vector_3;
typedef typename Tr::Bare_point Bare_point;
typedef typename Tr::Weighted_point Weighted_point;
public:
/**
* @brief Constructor
@ -445,11 +456,6 @@ protected:
std::vector<Vertex_handle>& modified_vertices,
bool *could_lock_zone = NULL) const
{
typedef typename C3T3::Triangulation Tr;
typedef typename Tr::Geom_traits Gt;
typedef typename Gt::FT FT;
typedef typename Tr::Bare_point Bare_point;
typedef Triangulation_helpers<typename C3T3::Triangulation> Th;
const Tr& tr = c3t3.triangulation();
@ -458,9 +464,9 @@ protected:
tr.geom_traits().compute_squared_length_3_object();
typename Gt::Construct_translated_point_3 translate =
tr.geom_traits().construct_translated_point_3_object();
typename Gt::Construct_point_3 wp2p =
typename Gt::Construct_point_3 cp =
tr.geom_traits().construct_point_3_object();
typename Gt::Construct_weighted_point_3 p2wp =
typename Gt::Construct_weighted_point_3 cwp =
tr.geom_traits().construct_weighted_point_3_object();
typename Gt::Construct_vector_3 vector =
tr.geom_traits().construct_vector_3_object();
@ -471,12 +477,12 @@ protected:
modified_vertices.clear();
// <periodic> shenanigans with v->point() ?
// norm depends on the local size of the mesh
FT sq_norm = this->compute_perturbation_sq_amplitude(v, c3t3, sq_step_size_);
FT step_length = CGAL::sqrt(sq_norm / sq_length(gradient_vector));
Vector_3 step_vector = step_length * gradient_vector;
Bare_point initial_loc = wp2p(v->point());
const Weighted_point& weighted_initial_loc = c3t3.triangulation().point(v);
Bare_point initial_loc = cp(weighted_initial_loc);
Bare_point new_loc = translate(initial_loc, step_vector);
Bare_point final_loc = new_loc;
@ -489,7 +495,7 @@ protected:
{
// as long as no topological change takes place
while(Th().no_topological_change__without_set_point(c3t3.triangulation(),
v, p2wp(final_loc))
v, cwp(final_loc))
&& (++i <= max_step_nb_) )
{
new_loc = translate(new_loc, step_vector);
@ -504,7 +510,7 @@ protected:
{
Vector_3 move_vector = vector(initial_loc, final_loc);
while( Th().no_topological_change(c3t3.triangulation(), v,
move_vector, p2wp(final_loc)) &&
move_vector, cwp(final_loc)) &&
++i <= max_step_nb_ )
{
new_loc = translate(new_loc, step_vector);
@ -525,7 +531,7 @@ protected:
// we know that there will be a combinatorial change
return helper.update_mesh_topo_change(v,
p2wp(final_loc),
cwp(final_loc),
criterion,
std::back_inserter(modified_vertices),
could_lock_zone);
@ -546,16 +552,21 @@ class Sq_radius_perturbation
: public Gradient_based_perturbation<C3T3,MeshDomain,SliverCriterion>
{
protected:
typedef Gradient_based_perturbation<C3T3,MeshDomain,SliverCriterion> Base;
typedef typename Base::Vertex_handle Vertex_handle;
typedef typename Base::Cell_handle Cell_handle;
typedef typename Base::FT FT;
typedef typename C3T3::Triangulation Tr;
typedef typename Tr::Geom_traits Gt;
typedef typename Gt::Vector_3 Vector_3;
typedef typename Tr::Bare_point Bare_point;
typedef typename Gt::Construct_point_3 Construct_point_3;
typedef Gradient_based_perturbation<C3T3, MeshDomain, SliverCriterion> Base;
typedef typename C3T3::Triangulation Tr;
typedef typename C3T3::Triangulation Triangulation;
typedef typename Base::Vertex_handle Vertex_handle;
typedef typename Base::Cell_handle Cell_handle;
typedef typename Tr::Geom_traits Gt;
typedef typename Gt::FT FT;
typedef typename Gt::Vector_3 Vector_3;
typedef typename Tr::Bare_point Bare_point;
typedef typename Tr::Weighted_point Weighted_point;
public:
/**
* @brief Constructor
@ -649,22 +660,25 @@ private:
const Cell_handle& cell,
const Vertex_handle& v) const
{
typedef typename C3T3::Triangulation Triangulation;
typedef typename Triangulation::Geom_traits Gt;
const Triangulation& tr = c3t3.triangulation();
typename Gt::Construct_translated_point_3 translate =
c3t3.triangulation().geom_traits().construct_translated_point_3_object();
typename Gt::Construct_point_3 wp2p =
typename Gt::Construct_point_3 cp =
c3t3.triangulation().geom_traits().construct_point_3_object();
// translate the tet so that cell->vertex((i+3)&3) is 0_{R^3}
unsigned int index = cell->index(v);
Vector_3 translate_to_origin(CGAL::ORIGIN, wp2p(tr.point(cell, (index+3)&3))); //p4
const Bare_point& p1 = translate(wp2p(tr.point(cell, index)), - translate_to_origin);
const Bare_point& p2 = translate(wp2p(tr.point(cell, (index+1)&3)), - translate_to_origin);
const Bare_point& p3 = translate(wp2p(tr.point(cell, (index+2)&3)), - translate_to_origin);
const Weighted_point& wvp = tr.point(cell, index);
const Weighted_point& wp2 = tr.point(cell, (index+1)&3);
const Weighted_point& wp3 = tr.point(cell, (index+2)&3);
const Weighted_point& wp4 = tr.point(cell, (index+3)&3);
// translate the tet so that 'wp4' is the origin
Vector_3 translate_to_origin(CGAL::ORIGIN, cp(wp4));
const Bare_point& p1 = translate(cp(wvp), - translate_to_origin);
const Bare_point& p2 = translate(cp(wp2), - translate_to_origin);
const Bare_point& p3 = translate(cp(wp3), - translate_to_origin);
// pre-compute everything
FT sq_p1 = p1.x()*p1.x() + p1.y()*p1.y() + p1.z()*p1.z();
@ -717,15 +731,19 @@ class Volume_perturbation
: public Gradient_based_perturbation<C3T3,MeshDomain,SliverCriterion>
{
protected:
typedef Gradient_based_perturbation<C3T3,MeshDomain,SliverCriterion> Base;
typedef typename Base::Vertex_handle Vertex_handle;
typedef typename Base::Cell_handle Cell_handle;
typedef typename Base::FT FT;
typedef typename C3T3::Triangulation::Geom_traits Gt;
typedef typename Gt::Vector_3 Vector_3;
typedef typename C3T3::Triangulation::Bare_point Bare_point;
typedef Gradient_based_perturbation<C3T3, MeshDomain, SliverCriterion> Base;
typedef typename Base::Vertex_handle Vertex_handle;
typedef typename Base::Cell_handle Cell_handle;
typedef typename C3T3::Triangulation Tr;
typedef typename Tr::Geom_traits Gt;
typedef typename Gt::FT FT;
typedef typename Gt::Vector_3 Vector_3;
typedef typename Tr::Bare_point Bare_point;
typedef typename Tr::Weighted_point Weighted_point;
public:
/**
* @brief Constructor
@ -823,7 +841,7 @@ private:
CGAL_assertion(cell->has_vertex(v));
const typename C3T3::Triangulation& tr = c3t3.triangulation();
typename Gt::Construct_point_3 wp2p =
typename Gt::Construct_point_3 cp =
c3t3.triangulation().geom_traits().construct_point_3_object();
const int i = cell->index(v);
@ -834,11 +852,14 @@ private:
if ( (i&1) == 0 )
std::swap(k1,k3);
const Bare_point& p1 = wp2p(tr.point(cell, k1));
const Bare_point& p2 = wp2p(tr.point(cell, k2));
const Bare_point& p3 = wp2p(tr.point(cell, k3));
const Weighted_point& wp1 = tr.point(cell, k1);
const Weighted_point& wp2 = tr.point(cell, k2);
const Weighted_point& wp3 = tr.point(cell, k3);
const Bare_point& p1 = cp(wp1);
const Bare_point& p2 = cp(wp2);
const Bare_point& p3 = cp(wp3);
FT gx = p2.y()*p3.z() + p1.y()*(p2.z()-p3.z())
- p3.y()*p2.z() - p1.z()*(p2.y()-p3.y());
@ -863,15 +884,18 @@ class Dihedral_angle_perturbation
: public Gradient_based_perturbation<C3T3,MeshDomain,SliverCriterion>
{
protected:
typedef Gradient_based_perturbation<C3T3,MeshDomain,SliverCriterion> Base;
typedef typename Base::Vertex_handle Vertex_handle;
typedef typename Base::Cell_handle Cell_handle;
typedef typename Base::FT FT;
typedef typename C3T3::Triangulation::Geom_traits Gt;
typedef typename Gt::Vector_3 Vector_3;
typedef typename C3T3::Triangulation::Bare_point Bare_point;
typedef typename Gt::Construct_point_3 Construct_point_3;
typedef Gradient_based_perturbation<C3T3, MeshDomain, SliverCriterion> Base;
typedef typename Base::Vertex_handle Vertex_handle;
typedef typename Base::Cell_handle Cell_handle;
typedef typename C3T3::Triangulation Tr;
typedef typename Tr::Geom_traits Gt;
typedef typename Gt::FT FT;
typedef typename Gt::Vector_3 Vector_3;
typedef typename Tr::Bare_point Bare_point;
typedef typename Tr::Weighted_point Weighted_point;
public:
/**
@ -970,13 +994,14 @@ private:
CGAL_assertion(cell->has_vertex(v));
const typename C3T3::Triangulation& tr = c3t3.triangulation();
typename Gt::Construct_point_3 wp2p =
typename Gt::Construct_point_3 cp =
tr.geom_traits().construct_point_3_object();
typename Gt::Compute_squared_distance_3 sq_distance =
tr.geom_traits().compute_squared_distance_3_object();
const int i = cell->index(v);
const Bare_point& p0 = wp2p(tr.point(cell, i));
const Weighted_point& wp0 = tr.point(cell, i);
const Bare_point& p0 = cp(wp0);
// Other indices
int k1 = (i+1)&3;
@ -995,9 +1020,12 @@ private:
std::swap(k2,k3);
// Here edge k1k2 minimizes dihedral angle
const Bare_point& p1 = wp2p(tr.point(cell, k1));
const Bare_point& p2 = wp2p(tr.point(cell, k2));
const Bare_point& p3 = wp2p(tr.point(cell, k3));
const Weighted_point& wp1 = tr.point(cell, k1);
const Weighted_point& wp2 = tr.point(cell, k2);
const Weighted_point& wp3 = tr.point(cell, k3);
const Bare_point& p1 = cp(wp1);
const Bare_point& p2 = cp(wp2);
const Bare_point& p3 = cp(wp3);
// grad of min dihedral angle (in cell) wrt p0
const Vector_3 p1p0 (p1,p0);
@ -1030,15 +1058,16 @@ private:
typename Gt::Compute_approximate_dihedral_angle_3 approx_dihedral_angle =
tr.geom_traits().compute_approximate_dihedral_angle_3_object();
typename Gt::Construct_point_3 wp2p =
typename Gt::Construct_point_3 cp =
tr.geom_traits().construct_point_3_object();
return CGAL::abs(approx_dihedral_angle(wp2p(p),
wp2p(tr.point(cell, k1)),
wp2p(tr.point(cell, k2)),
wp2p(tr.point(cell, k3))));
const Weighted_point& wp1 = tr.point(cell, k1);
const Weighted_point& wp2 = tr.point(cell, k2);
const Weighted_point& wp3 = tr.point(cell, k3);
return CGAL::abs(approx_dihedral_angle(cp(p), cp(wp1), cp(wp2), cp(wp3)));
}
/**
* @brief returns the cotangent of \c value
*/
@ -1055,7 +1084,7 @@ private:
{
const typename C3T3::Triangulation& tr = c3t3.triangulation();
typename Gt::Construct_point_3 wp2p = tr.geom_traits().construct_point_3_object();
typename Gt::Construct_point_3 cp = tr.geom_traits().construct_point_3_object();
int k1 = (i+1)&3;
int k2 = (i+2)&3;
@ -1065,9 +1094,12 @@ private:
if ( (i&1) == 1 )
std::swap(k1,k2);
const Bare_point& p1 = wp2p(tr.point(ch, k1));
const Bare_point& p2 = wp2p(tr.point(ch, k2));
const Bare_point& p3 = wp2p(tr.point(ch, k3));
const Weighted_point& wp1 = tr.point(ch, k1);
const Weighted_point& wp2 = tr.point(ch, k2);
const Weighted_point& wp3 = tr.point(ch, k3);
const Bare_point& p1 = cp(wp1);
const Bare_point& p2 = cp(wp2);
const Bare_point& p3 = cp(wp3);
// compute normal and return it
return normal(p1, p2, p3);
@ -1087,12 +1119,18 @@ class Random_based_perturbation
: public Abstract_perturbation<C3T3,MeshDomain,SliverCriterion>
{
protected:
typedef Abstract_perturbation<C3T3,MeshDomain,SliverCriterion> Base;
typedef typename Base::Vertex_handle Vertex_handle;
typedef typename Base::Cell_handle Cell_handle;
typedef typename Base::FT FT;
typedef typename C3T3::Triangulation::Geom_traits::Vector_3 Vector_3;
typedef Abstract_perturbation<C3T3, MeshDomain, SliverCriterion> Base;
typedef typename Base::Vertex_handle Vertex_handle;
typedef typename Base::Cell_handle Cell_handle;
typedef typename C3T3::Triangulation Tr;
typedef typename Tr::Geom_traits Gt;
typedef typename Gt::FT FT;
typedef typename Gt::Vector_3 Vector_3;
typedef typename Tr::Bare_point Bare_point;
public:
/**
* @brief constructor
@ -1150,8 +1188,6 @@ protected:
Vector_3 random_vector_fixed_size(const C3T3& c3t3,
const FT& vector_sq_size) const
{
typedef typename C3T3::Triangulation::Geom_traits Gt;
typename Gt::Compute_squared_length_3 sq_length =
c3t3.triangulation().geom_traits().compute_squared_length_3_object();
@ -1212,15 +1248,19 @@ class Li_random_perturbation
: public Random_based_perturbation<C3T3,MeshDomain,SliverCriterion>
{
protected:
typedef Random_based_perturbation<C3T3,MeshDomain,SliverCriterion> Base;
typedef typename Base::Vertex_handle Vertex_handle;
typedef typename Base::Cell_handle Cell_handle;
typedef typename Base::FT FT;
typedef typename C3T3::Triangulation::Geom_traits Gt;
typedef typename Gt::Vector_3 Vector_3;
typedef typename C3T3::Triangulation::Bare_point Bare_point;
typedef Random_based_perturbation<C3T3, MeshDomain, SliverCriterion> Base;
typedef typename Base::Vertex_handle Vertex_handle;
typedef typename Base::Cell_handle Cell_handle;
typedef typename C3T3::Triangulation Tr;
typedef typename Tr::Geom_traits Gt;
typedef typename Gt::FT FT;
typedef typename Gt::Vector_3 Vector_3;
typedef typename Tr::Bare_point Bare_point;
typedef typename Tr::Weighted_point Weighted_point;
public:
/**
* @brief constructor
@ -1288,16 +1328,13 @@ private:
{
typedef Triangulation_helpers<typename C3T3::Triangulation> Th;
typedef typename C3T3::Triangulation Tr;
typedef typename Tr::Geom_traits Gt;
typename Gt::Construct_translated_point_3 translate =
c3t3.triangulation().geom_traits().construct_translated_point_3_object();
typename Gt::Construct_weighted_point_3 p2wp =
typename Gt::Construct_weighted_point_3 cwp =
c3t3.triangulation().geom_traits().construct_weighted_point_3_object();
typename Gt::Construct_vector_3 vector =
c3t3.triangulation().geom_traits().construct_vector_3_object();
typename Gt::Construct_point_3 wp2p =
typename Gt::Construct_point_3 cp =
c3t3.triangulation().geom_traits().construct_point_3_object();
modified_vertices.clear();
@ -1308,8 +1345,9 @@ private:
// norm depends on the local size of the mesh
FT sq_norm = this->compute_perturbation_sq_amplitude(v, c3t3, this->sphere_sq_radius());
const Bare_point initial_location = wp2p(v->point());
const Weighted_point& weighted_initial_location = c3t3.triangulation().point(v);
const Bare_point initial_location = cp(weighted_initial_location);
// Initialize loop variables
bool criterion_improved = false;
Vertex_handle moving_vertex = v;
@ -1338,7 +1376,7 @@ private:
std::pair<bool,Vertex_handle> update =
helper.update_mesh(moving_vertex,
vector(initial_location, new_location),
p2wp(new_location),
cwp(new_location),
criterion,
std::back_inserter(tmp_mod_vertices),
could_lock_zone);
@ -1352,7 +1390,7 @@ private:
if ( update.first )
{
criterion_improved = true;
best_location = wp2p(moving_vertex->point());
best_location = cp(moving_vertex->point()); // <periodic> this seems useless
mod_vertices.insert(tmp_mod_vertices.begin(), tmp_mod_vertices.end());

View File

@ -618,10 +618,10 @@ Mesh_complex_3_in_triangulation_3(const Self& rhs)
const Vertex_handle& vb = it->left;
Vertex_handle new_va;
this->triangulation().is_vertex(va->point(), new_va);
this->triangulation().is_vertex(rhs.triangulation().point(va), new_va);
Vertex_handle new_vb;
this->triangulation().is_vertex(vb->point(), new_vb);
this->triangulation().is_vertex(rhs.triangulation().point(vb), new_vb);
this->add_to_complex(make_internal_edge(new_va,new_vb), it->info);
}
@ -631,7 +631,7 @@ Mesh_complex_3_in_triangulation_3(const Self& rhs)
end = rhs.corners_.end() ; it != end ; ++it )
{
Vertex_handle new_v;
this->triangulation().is_vertex(it->first->point(), new_v);
this->triangulation().is_vertex(rhs.triangulation().point(it->first), new_v);
this->add_to_complex(new_v, it->second);
}
@ -685,9 +685,10 @@ bool
Mesh_complex_3_in_triangulation_3<Tr,CI_,CSI_>::
is_valid(bool verbose) const
{
typedef typename Tr::Point::Point Bare_point;
typedef typename Tr::Point::Weight Weight;
typedef Weight FT;
typedef typename Tr::Bare_point Bare_point;
typedef typename Tr::Weighted_point Weighted_point;
typedef typename Tr::Weighted_point::Weight Weight;
typedef Weight FT;
std::map<Vertex_handle, int> vertex_map;
@ -712,7 +713,7 @@ is_valid(bool verbose) const
{
if(verbose)
std::cerr << "Validity error: vertex " << (void*)(&*vit->first)
<< " (" << vit->first->point() << ") "
<< " (" << this->triangulation().point(vit->first) << ") "
<< "is not a corner (dimension " << vit->first->in_dimension()
<< ") but has " << vit->second << " neighbor(s)!\n";
return false;
@ -727,11 +728,13 @@ is_valid(bool verbose) const
this->triangulation().geom_traits().construct_sphere_3_object();
typename Tr::Geom_traits::Do_intersect_3 do_intersect =
this->triangulation().geom_traits().do_intersect_3_object();
typename Tr::Geom_traits::Construct_point_3 wp2p =
typename Tr::Geom_traits::Construct_point_3 cp =
this->triangulation().geom_traits().construct_point_3_object();
const Bare_point& p = wp2p(it->right->point());
const Bare_point& q = wp2p(it->left->point());
const Weighted_point& itrwp = this->triangulation().point(it->right);
const Weighted_point& itlwp = this->triangulation().point(it->left);
const Bare_point& p = cp(itrwp);
const Bare_point& q = cp(itlwp);
const FT& sq_rp = it->right->point().weight();
const FT& sq_rq = it->left->point().weight();

View File

@ -512,9 +512,13 @@ public:
std::cout.flush();
#endif
CGAL::Random random(0);
typedef typename C3t3::Triangulation Tr;
typedef typename C3t3::Triangulation Tr;
typedef typename IGT::Sphere_3 Sphere_3;
typedef typename Polyhedron::Vertex_const_handle Vertex_const_handle;
Tr& tr = c3t3.triangulation();
typedef typename Polyhedron::Vertex_const_handle Vertex_const_handle;
typename Tr::Geom_traits::Construct_weighted_point_3 cwp
= tr.geom_traits().construct_weighted_point_3_object();
@ -538,10 +542,11 @@ public:
const Patch_id patch_id = vit->halfedge()->face()->patch_id();
CGAL_assertion(std::size_t(patch_id) <= nb_of_patch_plus_one);
typename Tr::Vertex_handle tr_v = tr.nearest_power_vertex(vit->point());
if (tr_v != typename Tr::Vertex_handle()) {
typedef typename IGT::Sphere_3 Sphere_3;
if (tr_v != typename Tr::Vertex_handle())
{
const Sphere_3 sphere(tr_v->point().point(), tr_v->point().weight());
if (!sphere.has_on_unbounded_side(vit->point())) continue;
if (!sphere.has_on_unbounded_side(vit->point()))
continue;
}
++nb_of_free_vertices_on_patch[patch_id];
}

View File

@ -64,7 +64,7 @@ namespace CGAL {
void operator()(const Vertex& vertex) const
{
typename Geom_traits::Construct_point_3 wp2p =
typename Geom_traits::Construct_point_3 cp =
r_c3t3_.triangulation().geom_traits().construct_point_3_object();
typename Geom_traits::Compute_weight_3 cw =
r_c3t3_.triangulation().geom_traits().compute_weight_3_object();
@ -72,7 +72,7 @@ namespace CGAL {
// Get vh properties
int dimension = vertex.in_dimension();
Weight w = (dimension < 2) ? cw(vertex.point()) : 0;
Weighted_point point(wp2p(vertex.point()), w);
Weighted_point point(cp(vertex.point()), w);
Index index = vertex.index();
// Insert point and restore handle properties

View File

@ -1,4 +1,8 @@
#define CGAL_MESH_3_VERBOSE
#define _DEBUG
#define CGAL_MESHES_DEBUG_REFINEMENT_POINTS
#define CGAL_MESH_3_DEBUG_FACET_CRITERIA
#define CGAL_MESH_3_DEBUG_CELL_CRITERIA
#include <CGAL/Mesh_3/config.h>
#include <CGAL/Periodic_3_mesh_3/config.h>
@ -72,8 +76,8 @@ int main(int argc, char** argv)
C3t3 c3t3 = CGAL::make_periodic_3_mesh_3<C3t3>(domain, criteria,
// odt(),
no_odt(),
// lloyd(),
no_lloyd(),
lloyd(),
// no_lloyd(),
// exude(),
no_exude(),
// perturb(sliver_bound=10, time_limit=10)

View File

@ -245,6 +245,11 @@ public:
Bare_point canonicalize_point(const Bare_point& p) const
{
if(p.x() >= domain().xmin() && p.x() < domain().xmax() &&
p.y() >= domain().ymin() && p.y() < domain().ymax() &&
p.z() >= domain().zmin() && p.z() < domain().zmax())
return p;
return robust_canonicalize_point(p);
}
@ -267,6 +272,11 @@ public:
Weighted_point canonicalize_point(const Weighted_point& p) const
{
if(p.x() >= domain().xmin() && p.x() < domain().xmax() &&
p.y() >= domain().ymin() && p.y() < domain().ymax() &&
p.z() >= domain().zmin() && p.z() < domain().zmax())
return p;
return robust_canonicalize_point(p);
}
@ -381,8 +391,13 @@ public:
// between p and q, for all combinations of offsets
FT min_squared_distance(const Bare_point& p, const Bare_point& q) const
{
// std::cout << "minsqd: " << p << " // " << q << std::endl;
const Bare_point cp = canonicalize_point(p);
const Bare_point cq = canonicalize_point(q);
// std::cout << "canon: " << cp << " // " << cq << std::endl;
FT min_sq_dist = std::numeric_limits<FT>::infinity();
for(int i = 0; i < 3; ++i) {
@ -397,6 +412,7 @@ public:
}
}
// std::cout << "minsqdt: " << min_sq_dist << std::endl;
return min_sq_dist;
}