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

Merge pull request #4976 from janetournois/Tet_remeshing-speedup-GF 

Tetrahedral Remeshing - speedup
This commit is contained in:
Sebastien Loriot 2020-10-02 15:29:24 +02:00 committed by GitHub
parent 47d2189f8e 8063d3cf6f
commit c54622e9b4
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
5 changed files with 387 additions and 93 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h
View File

@ -3648,10 +3648,17 @@ incident_slivers(const Vertex_handle& v,
std::vector<Cell_handle> incident_cells_;
tr_.incident_cells(v, std::back_inserter(incident_cells_));
#ifdef CGAL_CXX17
std::remove_copy_if(incident_cells_.begin(),
incident_cells_.end(),
out,
std::not_fn(Is_sliver<Sc>(c3t3_, criterion, sliver_bound)));
#else
std::remove_copy_if(incident_cells_.begin(),
incident_cells_.end(),
out,
std::not1(Is_sliver<Sc>(c3t3_,criterion,sliver_bound)));
#endif
return out;
}

235
Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/flip_edges.h
View File

@ -23,6 +23,7 @@
#include <boost/unordered_map.hpp>
#include <boost/unordered_set.hpp>
#include <boost/container/small_vector.hpp>
#include <boost/optional.hpp>
#include <limits>
#include <queue>
@ -86,11 +87,12 @@ void update_c3t3_facets(C3t3& c3t3,
}
}
template<typename C3t3>
template<typename C3t3, typename IncCellsVectorMap>
Sliver_removal_result flip_3_to_2(typename C3t3::Edge& edge,
C3t3& c3t3,
const std::vector<typename C3t3::Vertex_handle>& vertices_around_edge,
const Flip_Criterion& criterion)
const Flip_Criterion& criterion,
IncCellsVectorMap& inc_cells)
{
typedef typename C3t3::Triangulation Tr;
typedef typename C3t3::Facet Facet;
@ -186,19 +188,19 @@ Sliver_removal_result flip_3_to_2(typename C3t3::Edge& edge,
if (criterion == MIN_ANGLE_BASED)
{
//Current worst dihedral angle
FT curr_min_dh = min_dihedral_angle(tr, ch0);
curr_min_dh = (std::min)(curr_min_dh, min_dihedral_angle(tr, ch1));
curr_min_dh = (std::min)(curr_min_dh, min_dihedral_angle(tr, cell_to_remove));
Dihedral_angle_cosine curr_max_cosdh = max_cos_dihedral_angle(tr, ch0);
curr_max_cosdh = (std::max)(curr_max_cosdh, max_cos_dihedral_angle(tr, ch1));
curr_max_cosdh = (std::max)(curr_max_cosdh, max_cos_dihedral_angle(tr, cell_to_remove));
//Result worst dihedral angle
if (curr_min_dh > min_dihedral_angle(tr, vh2,
if (curr_max_cosdh < max_cos_dihedral_angle(tr, vh2,
ch0->vertex(indices(vh0_id, 0)),
ch0->vertex(indices(vh0_id, 1)),
ch0->vertex(indices(vh0_id, 2)))
|| curr_min_dh > min_dihedral_angle(tr, vh3,
ch1->vertex(indices(vh1_id, 0)),
ch1->vertex(indices(vh1_id, 1)),
ch1->vertex(indices(vh1_id, 2))))
|| curr_max_cosdh < max_cos_dihedral_angle(tr, vh3,
ch1->vertex(indices(vh1_id, 0)),
ch1->vertex(indices(vh1_id, 1)),
ch1->vertex(indices(vh1_id, 2))))
return NO_BEST_CONFIGURATION;
}
else if (criterion == AVERAGE_ANGLE_BASED)
@ -214,7 +216,7 @@ Sliver_removal_result flip_3_to_2(typename C3t3::Edge& edge,
(min_dihedral_angle(tr, vh2, ch0->vertex(indices(vh0_id, 0)),
ch0->vertex(indices(vh0_id, 1)),
ch0->vertex(indices(vh0_id, 2)))
+ min_dihedral_angle(tr, vh3, ch1->vertex(indices(vh1_id, 0)),
+ min_dihedral_angle(tr, vh3, ch1->vertex(indices(vh1_id, 0)),
ch1->vertex(indices(vh1_id, 1)),
ch1->vertex(indices(vh1_id, 2))));
//Result worst dihedral angle
@ -312,6 +314,9 @@ Sliver_removal_result flip_3_to_2(typename C3t3::Edge& edge,
ch->set_neighbor(v, mirror_facet.first);
}
ch->vertex(v)->set_cell(ch);
inc_cells[ch->vertex(v)] = boost::none;
ch->reset_cache_validity();
}
}
@ -372,7 +377,7 @@ void find_best_flip_to_improve_dh(C3t3& c3t3,
typename C3t3::Vertex_handle vh2,
typename C3t3::Vertex_handle vh3,
CandidatesQueue& candidates,
double curr_min_dh,
const Dihedral_angle_cosine& curr_max_cos_dh,
bool is_sliver_well_oriented = true,
int e_id = 0)
{
@ -382,8 +387,6 @@ void find_best_flip_to_improve_dh(C3t3& c3t3,
typedef typename C3t3::Facet Facet;
typedef typename Tr::Facet_circulator Facet_circulator;
typedef typename Tr::Cell_circulator Cell_circulator;
typedef typename Tr::Geom_traits Gt;
typedef typename Gt::FT FT;
// std::cout << "find_best_flip_to_improve_dh boundary " << std::endl;
Tr& tr = c3t3.triangulation();
@ -451,12 +454,12 @@ void find_best_flip_to_improve_dh(C3t3& c3t3,
Cell_circulator cell_circulator = tr.incident_cells(edge);
Cell_circulator done = cell_circulator;
boost::container::small_vector<Facet, 60> facets;
for (std::size_t i = 0; i < opposite_vertices.size(); ++i)
{
Vertex_handle vh = opposite_vertices[i];
bool keep = true;
std::vector<Facet> facets;
do
{
//Store it if it do not have vh
@ -474,7 +477,7 @@ void find_best_flip_to_improve_dh(C3t3& c3t3,
} while (++cell_circulator != done);
FT min_flip_dihedral_angle = (std::numeric_limits<FT>::max)();
Dihedral_angle_cosine max_flip_cos_dh(CGAL::NEGATIVE, 1., 1.);
for (const Facet& fi : facets)
{
if (!tr.is_infinite(fi.first))
@ -483,32 +486,75 @@ void find_best_flip_to_improve_dh(C3t3& c3t3,
fi.first->vertex(indices(fi.second, 1)),
fi.first->vertex(indices(fi.second, 2))))
{
min_flip_dihedral_angle = (std::min)(min_flip_dihedral_angle,
min_dihedral_angle(tr, vh, fi.first->vertex(indices(fi.second, 0)),
fi.first->vertex(indices(fi.second, 1)),
fi.first->vertex(indices(fi.second, 2))));
max_flip_cos_dh = (std::max)(
max_flip_cos_dh,
max_cos_dihedral_angle(tr, vh, fi.first->vertex(indices(fi.second, 0)),
fi.first->vertex(indices(fi.second, 1)),
fi.first->vertex(indices(fi.second, 2))));
}
else
{
keep = false;
break;
}
if (max_flip_cos_dh.is_one())//it will not get worse than 1.
{
keep = false;
break;
}
}
}
facets.clear();
if (keep && (curr_min_dh < min_flip_dihedral_angle || !is_sliver_well_oriented))
if (keep && (max_flip_cos_dh < curr_max_cos_dh || !is_sliver_well_oriented))
{
//std::cout << "vh " << vh->info() <<" old " << curr_min_dh << " min " << min_flip_dihedral_angle << std::endl;
candidates.push(std::make_pair(min_flip_dihedral_angle, std::make_pair(vh, e_id)));
//std::cout << "vh " << vh->info() <<" old " << curr_max_cos_dh << " min " << min_flip_tan_dh << std::endl;
candidates.push(std::make_pair(max_flip_cos_dh, std::make_pair(vh, e_id)));
}
}
}
template<typename C3t3, typename CandidatesQueue>
template<typename Vertex_handle, typename CellVector, typename Cell_handle>
bool is_edge_uv(Vertex_handle u,
Vertex_handle v,
const CellVector& cells_incident_to_u,
Cell_handle& cell,
int& i,
int& j)
{
if (u == v)
return false;
for (typename CellVector::value_type c : cells_incident_to_u)
{
if (c->has_vertex(v, j))
{
cell = c;
i = cell->index(u);
return true;
}
}
return false;
}
template<typename Vertex_handle, typename CellVector>
bool is_edge_uv(Vertex_handle u,
Vertex_handle v,
const CellVector& cells_incident_to_u)
{
typename CellVector::value_type c;
int i, j;
return is_edge_uv(u, v, cells_incident_to_u, c, i, j);
}
template<typename C3t3, typename CandidatesQueue,
typename IncCellsVectorMap>
void find_best_flip_to_improve_dh(C3t3& c3t3,
typename C3t3::Edge& edge,
CandidatesQueue& candidates,
double curr_min_dh,
const Dihedral_angle_cosine& curr_max_cosdh,
IncCellsVectorMap& inc_cells,
bool is_sliver_well_oriented = true,
int e_id = 0)
{
@ -518,8 +564,6 @@ void find_best_flip_to_improve_dh(C3t3& c3t3,
typedef typename C3t3::Facet Facet;
typedef typename Tr::Facet_circulator Facet_circulator;
typedef typename Tr::Cell_circulator Cell_circulator;
typedef typename Tr::Geom_traits Gt;
typedef typename Gt::FT FT;
Tr& tr = c3t3.triangulation();
@ -547,6 +591,17 @@ void find_best_flip_to_improve_dh(C3t3& c3t3,
}
}
if(tr.is_infinite(vh))
continue;
boost::optional<boost::container::small_vector<Cell_handle, 64>>& o_inc_vh = inc_cells[vh];
if (o_inc_vh == boost::none)
{
boost::container::small_vector<Cell_handle, 64> inc_vec;
tr.incident_cells(vh, std::back_inserter(inc_vec));
o_inc_vh = inc_vec;
}
Facet_circulator facet_circulator = curr_fcirc;
Facet_circulator facet_done = curr_fcirc;
@ -563,15 +618,16 @@ void find_best_flip_to_improve_dh(C3t3& c3t3,
indices(facet_circulator->second, i));
if (curr_vertex != vh0 && curr_vertex != vh1)
{
Cell_handle ch;
int i0, i1;
if (tr.is_edge(curr_vertex, vh, ch, i0, i1))
if (is_edge_uv(vh, curr_vertex, boost::get(o_inc_vh)))
{
is_edge = true;
break;
}
}
}
} while (++facet_circulator != facet_done);
if (!is_edge && !tr.is_infinite(vh))
if (!is_edge)
opposite_vertices.push_back(vh);
nb_cells_around_edge++;
@ -589,12 +645,10 @@ void find_best_flip_to_improve_dh(C3t3& c3t3,
Cell_circulator cell_circulator = tr.incident_cells(edge);
Cell_circulator done = cell_circulator;
for (std::size_t i = 0; i < opposite_vertices.size(); ++i)
boost::container::small_vector<Facet, 60> facets;
for (Vertex_handle vh : opposite_vertices)
{
Vertex_handle vh = opposite_vertices[i];
bool keep = true;
std::vector<Facet> facets;
do
{
//Store it if it do not have vh
@ -612,7 +666,7 @@ void find_best_flip_to_improve_dh(C3t3& c3t3,
}
while (++cell_circulator != done);
FT min_flip_dihedral_angle = (std::numeric_limits<FT>::max)();
Dihedral_angle_cosine max_flip_cos_dh(CGAL::NEGATIVE, 1., 1.);
for (const Facet& fi : facets)
{
if (!tr.is_infinite(fi.first))
@ -621,31 +675,39 @@ void find_best_flip_to_improve_dh(C3t3& c3t3,
fi.first->vertex(indices(fi.second, 1)),
fi.first->vertex(indices(fi.second, 2))))
{
min_flip_dihedral_angle = (std::min)(min_flip_dihedral_angle,
min_dihedral_angle(tr, vh, fi.first->vertex(indices(fi.second, 0)),
fi.first->vertex(indices(fi.second, 1)),
fi.first->vertex(indices(fi.second, 2))));
max_flip_cos_dh = (std::max)(max_flip_cos_dh,
max_cos_dihedral_angle(tr, vh, fi.first->vertex(indices(fi.second, 0)),
fi.first->vertex(indices(fi.second, 1)),
fi.first->vertex(indices(fi.second, 2))));
}
else
{
keep = false;
break;
}
if (max_flip_cos_dh.is_one())//it will not get worse than 1.
{
keep = false;
break;
}
}
}
facets.clear();
if (keep && (curr_min_dh < min_flip_dihedral_angle || !is_sliver_well_oriented))
if (keep && (max_flip_cos_dh < curr_max_cosdh || !is_sliver_well_oriented))
{
//std::cout << "vh " << vh->info() <<" old " << curr_min_dh << " min " << min_flip_dihedral_angle << std::endl;
candidates.push(std::make_pair(min_flip_dihedral_angle, std::make_pair(vh, e_id)));
//std::cout << "vh " << vh->info() <<" old " << curr_max_cosdh << " min " << min_flip_tan_dh << std::endl;
candidates.push(std::make_pair(max_flip_cos_dh, std::make_pair(vh, e_id)));
}
}
}
template<typename C3t3, typename Visitor>
template<typename C3t3, typename IncCellsVectorMap, typename Visitor>
Sliver_removal_result flip_n_to_m(C3t3& c3t3,
typename C3t3::Edge& edge,
typename C3t3::Vertex_handle vh,
IncCellsVectorMap& inc_cells,
Visitor& visitor,
bool check_validity = false)
{
@ -693,7 +755,14 @@ Sliver_removal_result flip_n_to_m(C3t3& c3t3,
facet_circulator++;
facet_circulator++;
std::vector<Vertex_handle> vertices_around_edge;
boost::optional<boost::container::small_vector<Cell_handle, 64>>& o_inc_vh = inc_cells[vh];
if (o_inc_vh == boost::none)
{
boost::container::small_vector<Cell_handle, 64> inc_vec;
tr.incident_cells(vh, std::back_inserter(inc_vec));
o_inc_vh = inc_vec;
}
do
{
//Get the ids of the opposite vertices
@ -703,18 +772,13 @@ Sliver_removal_result flip_n_to_m(C3t3& c3t3,
indices(facet_circulator->second, i));
if (curr_vertex != vh0 && curr_vertex != vh1)
{
Cell_handle ch;
int i0, i1;
if (tr.is_edge(curr_vertex, vh, ch, i0, i1))
if (is_edge_uv(vh, curr_vertex, boost::get(o_inc_vh)))
return NOT_FLIPPABLE;
vertices_around_edge.push_back(curr_vertex);
}
}
} while (++facet_circulator != facet_done);
std::vector<Cell_handle> cells_around_edge;
boost::container::small_vector<Cell_handle, 20> to_remove;
//Neighbors that will need to be updated after flip
@ -733,8 +797,6 @@ Sliver_removal_result flip_n_to_m(C3t3& c3t3,
Cell_circulator done = cell_circulator;
do
{
cells_around_edge.push_back(cell_circulator);
//Facets opposite to vh0
Facet facet_vh0(cell_circulator, cell_circulator->index(vh0));
neighbor_facets.insert(tr.mirror_facet(facet_vh0));
@ -877,6 +939,9 @@ Sliver_removal_result flip_n_to_m(C3t3& c3t3,
ch->set_neighbor(v, facet.first);
}
ch->vertex(v)->set_cell(ch);
inc_cells[ch->vertex(v)] = boost::none;
ch->reset_cache_validity();
}
}
@ -943,54 +1008,61 @@ Sliver_removal_result flip_n_to_m(C3t3& c3t3,
}
template<typename C3t3, typename Visitor>
template<typename C3t3, typename IncCellsVectorMap, typename Visitor>
Sliver_removal_result flip_n_to_m(typename C3t3::Edge& edge,
C3t3& c3t3,
std::vector<typename C3t3::Vertex_handle>& boundary_vertices,
const std::vector<typename C3t3::Vertex_handle>& boundary_vertices,
const Flip_Criterion& criterion,
IncCellsVectorMap& inc_cells,
Visitor& visitor)
{
typedef typename C3t3::Vertex_handle Vertex_handle;
typedef typename C3t3::Triangulation::Cell_circulator Cell_circulator;
typedef typename C3t3::Triangulation::Geom_traits Gt;
typedef typename Gt::FT FT;
typename C3t3::Triangulation& tr = c3t3.triangulation();
Sliver_removal_result result = NOT_FLIPPABLE;
typedef std::pair<double, std::pair<Vertex_handle, int> > Angle_and_vertex;
typedef std::pair<Dihedral_angle_cosine, std::pair<Vertex_handle, int> > CosAngle_and_vertex;
//std::cout << "n_to_m_flip " << boundary_vertices.size() << std::endl;
if (criterion == MIN_ANGLE_BASED)
{
std::priority_queue<Angle_and_vertex> candidates;
std::priority_queue<CosAngle_and_vertex,
std::vector<CosAngle_and_vertex>,
std::greater<CosAngle_and_vertex>
> candidates;
Cell_circulator circ = c3t3.triangulation().incident_cells(edge);
Cell_circulator done = circ;
FT curr_min_dh = min_dihedral_angle(tr, circ++);
Dihedral_angle_cosine curr_max_cosdh = max_cos_dihedral_angle(tr, circ++);
while (circ != done)
{
curr_min_dh = (std::min)(curr_min_dh, min_dihedral_angle(tr, circ++));
curr_max_cosdh = (std::max)(curr_max_cosdh, max_cos_dihedral_angle(tr, circ++));
}
if (boundary_vertices.size() == 2)
find_best_flip_to_improve_dh(c3t3, edge, boundary_vertices[0], boundary_vertices[1],
candidates, curr_min_dh);
candidates, curr_max_cosdh);
else
find_best_flip_to_improve_dh(c3t3, edge, candidates, curr_min_dh);
find_best_flip_to_improve_dh(c3t3, edge, candidates, curr_max_cosdh,
inc_cells);
bool flip_performed = false;
while (!candidates.empty() && !flip_performed)
{
Angle_and_vertex curr_cost_vpair = candidates.top();
CosAngle_and_vertex curr_cost_vpair = candidates.top();
candidates.pop();
//std::cout << curr_min_dh << " old, current " << curr_cost_vpair.second.first->info() <<" and angle " << curr_cost_vpair.first << std::endl;
// std::cout << "\tcurrent cos = " << curr_max_cosdh
// << "\t angle = " << std::acos(curr_max_cosdh) * 180./CGAL_PI << std::endl;
// std::cout << "\tcandidate cos = " << curr_cost_vpair.first
// << "\t angle = " << std::acos(curr_cost_vpair.first) * 180./CGAL_PI << std::endl;
// std::cout << std::endl;
if (curr_min_dh >= curr_cost_vpair.first)
if (curr_max_cosdh <= curr_cost_vpair.first)
return NO_BEST_CONFIGURATION;
result = flip_n_to_m(c3t3, edge, curr_cost_vpair.second.first, visitor);
result = flip_n_to_m(c3t3, edge, curr_cost_vpair.second.first, inc_cells, visitor);
if (result != NOT_FLIPPABLE)
flip_performed = true;
@ -1000,10 +1072,11 @@ Sliver_removal_result flip_n_to_m(typename C3t3::Edge& edge,
return result;
}
template<typename C3t3, typename Visitor>
template<typename C3t3, typename IncCellsVectorMap, typename Visitor>
Sliver_removal_result find_best_flip(typename C3t3::Edge& edge,
C3t3& c3t3,
const Flip_Criterion& criterion,
IncCellsVectorMap& inc_cells,
Visitor& visitor)
{
typedef typename C3t3::Triangulation Tr;
@ -1066,7 +1139,7 @@ Sliver_removal_result find_best_flip(typename C3t3::Edge& edge,
{
std::vector<Vertex_handle> vertices;
vertices.insert(vertices.end(), vertices_around_edge.begin(), vertices_around_edge.end());
res = flip_3_to_2(edge, c3t3, vertices, criterion);
res = flip_3_to_2(edge, c3t3, vertices, criterion, inc_cells);
}
}
else
@ -1078,12 +1151,11 @@ Sliver_removal_result find_best_flip(typename C3t3::Edge& edge,
{
std::vector<Vertex_handle> vertices;
vertices.insert(vertices.end(), boundary_vertices.begin(), boundary_vertices.end());
res = flip_n_to_m(edge, c3t3, vertices, criterion, visitor);
res = flip_n_to_m(edge, c3t3, vertices, criterion, inc_cells, visitor);
//return n_to_m_flip(edge, boundary_vertices, flip_criterion);
}
}
return res;
}
@ -1096,20 +1168,33 @@ std::size_t flip_all_edges(const std::vector<VertexPair>& edges,
{
typedef typename C3t3::Triangulation Tr;
typedef typename Tr::Cell_handle Cell_handle;
typedef typename Tr::Vertex_handle Vertex_handle;
typedef typename Tr::Edge Edge;
Tr& tr = c3t3.triangulation();
std::unordered_map<Vertex_handle,
boost::optional<boost::container::small_vector<Cell_handle, 64> > > inc_cells;
std::size_t count = 0;
for (const VertexPair vp : edges)
{
boost::optional<boost::container::small_vector<Cell_handle, 64>>&
o_inc_vh = inc_cells[vp.first];
if (o_inc_vh == boost::none)
{
boost::container::small_vector<Cell_handle, 64> inc_vec;
tr.incident_cells(vp.first, std::back_inserter(inc_vec));
o_inc_vh = inc_vec;
}
Cell_handle ch;
int i0, i1;
if (tr.is_edge(vp.first, vp.second, ch, i0, i1))
if (is_edge_uv(vp.first, vp.second, boost::get(o_inc_vh), ch, i0, i1))
{
Edge edge(ch, i0, i1);
Sliver_removal_result res = find_best_flip(edge, c3t3, criterion, visitor);
Sliver_removal_result res = find_best_flip(edge, c3t3, criterion, inc_cells, visitor);
if (res == INVALID_CELL || res == INVALID_VERTEX || res == INVALID_ORIENTATION)
{
std::cout << "FLIP PROBLEM!!!!" << std::endl;
@ -1126,6 +1211,10 @@ std::size_t flip_all_edges(const std::vector<VertexPair>& edges,
}
}
}
for(Cell_handle c : c3t3.triangulation().finite_cell_handles())
c->reset_cache_validity();
return count;
}

47
Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/smooth_vertices.h
View File

@ -328,7 +328,12 @@ private:
bool check_inversion_and_move(const typename Tr::Vertex_handle v,
const typename Tr::Point& final_pos,
const CellRange& inc_cells,
const Tr& /* tr */)
const Tr& /* tr */,
#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE
double& total_move)
#else
double&)
#endif
{
const typename Tr::Point backup = v->point(); //backup v's position
const typename Tr::Geom_traits::Point_3 pv = point(backup);
@ -336,6 +341,7 @@ private:
bool valid_orientation = false;
double frac = 1.0;
typename Tr::Geom_traits::Vector_3 move(pv, point(final_pos));
do
{
v->set_point(typename Tr::Point(pv + frac * move));
@ -354,15 +360,17 @@ private:
}
}
valid_orientation = valid_try;
// std::cout << std::boolalpha << "valid orientation = " << valid_orientation
// << "\tfrac = " << frac << std::endl;
}
while(!valid_orientation && frac > 0.1);
if (!valid_orientation) //move failed
v->set_point(backup);
#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE
else
total_move += CGAL::approximate_sqrt(CGAL::squared_distance(pv, point(v->point())));
#endif
return valid_orientation;
}
@ -407,7 +415,10 @@ public:
std::cout << "Smooth vertices...";
std::cout.flush();
#endif
std::size_t nb_done = 0;
std::size_t nb_done_3d = 0;
std::size_t nb_done_2d = 0;
std::size_t nb_done_1d = 0;
FT total_move = 0.;
Tr& tr = c3t3.triangulation();
@ -521,8 +532,8 @@ public:
#endif
// move vertex
const typename Tr::Point new_pos(final_position.x(), final_position.y(), final_position.z());
if(check_inversion_and_move(v, new_pos, inc_cells[vid], tr))
nb_done++;
if(check_inversion_and_move(v, new_pos, inc_cells[vid], tr, total_move))
nb_done_1d++;
}
else if (neighbors[vid] > 0)
{
@ -555,8 +566,8 @@ public:
#endif
// move vertex
const typename Tr::Point new_pos(final_position.x(), final_position.y(), final_position.z());
if(check_inversion_and_move(v, new_pos, inc_cells[vid], tr))
nb_done++;
if(check_inversion_and_move(v, new_pos, inc_cells[vid], tr, total_move))
nb_done_1d++;
}
}
}
@ -628,8 +639,8 @@ public:
os_surf << "2 " << current_pos << " " << final_position << std::endl;
#endif
const typename Tr::Point new_pos(final_position.x(), final_position.y(), final_position.z());
if(check_inversion_and_move(v, new_pos, inc_cells[vid], tr))
nb_done++;
if(check_inversion_and_move(v, new_pos, inc_cells[vid], tr, total_move))
nb_done_2d++;
}
else if (neighbors[vid] > 0)
{
@ -641,8 +652,8 @@ public:
if (boost::optional<Vector_3> mls_projection = project(si, current_pos))
{
const typename Tr::Point new_pos(CGAL::ORIGIN + *mls_projection);
if(check_inversion_and_move(v, new_pos, inc_cells[vid], tr))
nb_done++;
if(check_inversion_and_move(v, new_pos, inc_cells[vid], tr, total_move))
nb_done_2d++;
#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
os_surf0 << "2 " << current_pos << " " << new_pos << std::endl;
@ -693,8 +704,8 @@ public:
#endif
const Vector_3 p = smoothed_positions[vid] / static_cast<FT>(neighbors[vid]);
typename Tr::Point new_pos(p.x(), p.y(), p.z());
if(check_inversion_and_move(v, new_pos, inc_cells[vid], tr))
nb_done++;
if(check_inversion_and_move(v, new_pos, inc_cells[vid], tr, total_move))
nb_done_3d++;
#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
os_vol << " " << point(v->point()) << std::endl;
@ -704,7 +715,11 @@ public:
CGAL_assertion(CGAL::Tetrahedral_remeshing::debug::are_cell_orientations_valid(tr));
#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE
std::cout << " done (" << nb_done << " vertices smoothed)." << std::endl;
std::size_t nb_done = nb_done_3d + nb_done_2d + nb_done_1d;
std::cout << " done ("
<< nb_done_3d << "/" << nb_done_2d << "/" << nb_done_1d << " vertices smoothed,"
<< " average move = " << (total_move / nb_done)
<< ")." << std::endl;
#endif
#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
CGAL::Tetrahedral_remeshing::debug::dump_vertices_by_dimension(

14
Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/tetrahedral_adaptive_remeshing_impl.h
View File

@ -368,8 +368,18 @@ public:
CGAL::Tetrahedral_remeshing::debug::dump_c3t3(m_c3t3, "99-postprocess");
#endif
#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE
std::cout << "(peeling removed " << nb_slivers_peel << " slivers)" << std::endl;
std::cout << "done." << std::endl;
mindh = 180.;
for (Cell_handle cit : tr().finite_cell_handles())
{
if (m_c3t3.is_in_complex(cit))
{
const double dh = min_dihedral_angle(tr(), cit);
mindh = (std::min)(dh, mindh);
}
}
std::cout << "Peeling done (removed " << nb_slivers_peel << " slivers, "
<< "min dihedral angle = " << mindh << ")." << std::endl;
#endif
return nb_slivers_peel;
}

177
Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/tetrahedral_remeshing_helpers.h
View File

@ -98,7 +98,7 @@ typename Geom_traits::FT min_dihedral_angle(const Point& p,
FT a = CGAL::abs(dihedral_angle(p, q, r, s, gt));
FT min_dh = a;
a = CGAL::abs(dihedral_angle(p, r, q, s, gt));
a = CGAL::abs(dihedral_angle(p, r, s, q, gt));
min_dh = (std::min)(a, min_dh);
a = CGAL::abs(dihedral_angle(p, s, q, r, gt));
@ -107,7 +107,7 @@ typename Geom_traits::FT min_dihedral_angle(const Point& p,
a = CGAL::abs(dihedral_angle(q, r, p, s, gt));
min_dh = (std::min)(a, min_dh);
a = CGAL::abs(dihedral_angle(q, s, p, r, gt));
a = CGAL::abs(dihedral_angle(q, s, r, p, gt));
min_dh = (std::min)(a, min_dh);
a = CGAL::abs(dihedral_angle(r, s, p, q, gt));
@ -141,6 +141,179 @@ typename Tr::Geom_traits::FT min_dihedral_angle(const Tr& tr,
c->vertex(3));
}
struct Dihedral_angle_cosine
{
CGAL::Sign m_sgn;
double m_sq_num;
double m_sq_den;
Dihedral_angle_cosine(const CGAL::Sign& sgn, const double& sq_num, const double& sq_den)
: m_sgn(sgn)
, m_sq_num(sq_num)
, m_sq_den(sq_den)
{}
bool is_one() const
{
return m_sgn == CGAL::POSITIVE && m_sq_num == m_sq_den;
}
double signed_square_value() const
{
switch(m_sgn)
{
case CGAL::POSITIVE:
return m_sq_num / m_sq_den;
case CGAL::ZERO:
return 0.;
default:
CGAL_assertion(m_sgn == CGAL::NEGATIVE);
return -1. * m_sq_num / m_sq_den;
};
}
friend bool operator<(const Dihedral_angle_cosine& l,
const Dihedral_angle_cosine& r)
{
//if numerators have different signs
if (l.m_sgn == CGAL::NEGATIVE && r.m_sgn != CGAL::NEGATIVE)
return true;
else if (l.m_sgn == CGAL::POSITIVE && r.m_sgn != CGAL::POSITIVE)
return false;
else if (l.m_sgn == CGAL::ZERO)
return (r.m_sgn == CGAL::POSITIVE);
//else numerators have the same sign
else if (l.m_sgn == CGAL::POSITIVE) //both angles are in [0; PI/2[
{
CGAL_assertion(r.m_sgn == CGAL::POSITIVE);
return (l.m_sq_num * r.m_sq_den < r.m_sq_num* l.m_sq_den);
}
else //both angles are in [PI/2; PI]
{
CGAL_assertion(l.m_sgn != CGAL::POSITIVE);
CGAL_assertion(r.m_sgn != CGAL::POSITIVE);
return (l.m_sq_num * r.m_sq_den >= r.m_sq_num* l.m_sq_den);
}
}
friend bool operator<=(const Dihedral_angle_cosine& l,
const Dihedral_angle_cosine& r)
{
if(l < r)
return true;
else
return l.m_sgn == r.m_sgn
&& l.m_sq_num * r.m_sq_den == r.m_sq_num * l.m_sq_den;
}
};
template<typename Gt>
Dihedral_angle_cosine cos_dihedral_angle(const typename Gt::Point_3& i,
const typename Gt::Point_3& j,
const typename Gt::Point_3& k,
const typename Gt::Point_3& l,
const Gt& gt)
{
CGAL_assertion(CGAL::orientation(i, j, k, l) != CGAL::NEGATIVE);
typename Gt::Construct_vector_3 vector = gt.construct_vector_3_object();
typename Gt::Construct_cross_product_vector_3 cross_product =
gt.construct_cross_product_vector_3_object();
typename Gt::Compute_scalar_product_3 scalar_product =
gt.compute_scalar_product_3_object();
typedef typename Gt::FT FT;
typedef typename Gt::Vector_3 Vector_3;
const Vector_3 ij = vector(i, j);
const Vector_3 ik = vector(i, k);
const Vector_3 il = vector(i, l);
const Vector_3 ijik = cross_product(ij, ik);
if(CGAL::NULL_VECTOR == ijik)
return Dihedral_angle_cosine(CGAL::POSITIVE, 1.,1.);
const Vector_3 ilij = cross_product(il, ij);
if (CGAL::NULL_VECTOR == ilij)
return Dihedral_angle_cosine(CGAL::POSITIVE, 1.,1.);
const FT num = scalar_product(ijik, ilij);
if(num == 0.)
return Dihedral_angle_cosine(CGAL::ZERO, 0.,1.);
const double sqden = CGAL::to_double(
scalar_product(ijik, ijik) * scalar_product(ilij, ilij));
return Dihedral_angle_cosine(CGAL::sign(num), CGAL::square(num), sqden);
}
template<typename Point, typename Geom_traits>
Dihedral_angle_cosine max_cos_dihedral_angle(const Point& p,
const Point& q,
const Point& r,
const Point& s,
const Geom_traits& gt)
{
Dihedral_angle_cosine a = cos_dihedral_angle(p, q, r, s, gt);
Dihedral_angle_cosine max_cos_dh = a;
if(max_cos_dh.is_one()) return max_cos_dh;
a = cos_dihedral_angle(p, r, s, q, gt);
if(max_cos_dh < a) max_cos_dh = a;
if (max_cos_dh.is_one()) return max_cos_dh;
a = cos_dihedral_angle(p, s, q, r, gt);
if (max_cos_dh < a) max_cos_dh = a;
if (max_cos_dh.is_one()) return max_cos_dh;
a = cos_dihedral_angle(q, r, p, s, gt);
if (max_cos_dh < a) max_cos_dh = a;
if (max_cos_dh.is_one()) return max_cos_dh;
a = cos_dihedral_angle(q, s, r, p, gt);
if (max_cos_dh < a) max_cos_dh = a;
if (max_cos_dh.is_one()) return max_cos_dh;
a = cos_dihedral_angle(r, s, p, q, gt);
if (max_cos_dh < a) max_cos_dh = a;
return max_cos_dh;
}
template<typename Tr>
Dihedral_angle_cosine max_cos_dihedral_angle(const Tr& tr,
const typename Tr::Vertex_handle v0,
const typename Tr::Vertex_handle v1,
const typename Tr::Vertex_handle v2,
const typename Tr::Vertex_handle v3)
{
return max_cos_dihedral_angle(point(v0->point()),
point(v1->point()),
point(v2->point()),
point(v3->point()),
tr.geom_traits());
}
template<typename Tr>
Dihedral_angle_cosine max_cos_dihedral_angle(const Tr& tr,
const typename Tr::Cell_handle c)
{
if (c->is_cache_valid())
return Dihedral_angle_cosine(CGAL::sign(c->sliver_value()),
CGAL::abs(c->sliver_value()), 1.);
Dihedral_angle_cosine cos_dh = max_cos_dihedral_angle(tr,
c->vertex(0),
c->vertex(1),
c->vertex(2),
c->vertex(3));
c->set_sliver_value(cos_dh.signed_square_value());
return cos_dh;
}
template<typename C3t3>
bool is_peelable(const C3t3& c3t3,
const typename C3t3::Cell_handle ch,