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

Merge pull request #7535 from lrineau/Triangulation_3-fix_simplex_traverser-GF-CGAL-5.6 

fix simplex traverser for CGAL 5.6
This commit is contained in:
Laurent Rineau 2023-07-05 16:32:36 +02:00
parent 2e8ce4006d 7e2e444825
commit 04077e9f4f
8 changed files with 1116 additions and 715 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
STL_Extension/include/CGAL/Iterator_range.h
View File

@ -78,6 +78,13 @@ namespace CGAL {
{
return std::tuple<const I&, const I&>{this->first, this->second};
}
template <template<class...> class Container>
auto to() const
{
using V = std::remove_cv_t<std::remove_reference_t<decltype(*begin())>>;
return Container<V>(begin(), end());
}
};
template <typename T>

7
TDS_3/include/CGAL/Triangulation_simplex_3.h
View File

@ -96,8 +96,10 @@ public:
// returns the dimension of the simplex
int dimension () const {
return (ref & 3);
if(ref == -1) return -1;
else return (ref & 3);
}
// returns an incident cell:
Cell_handle incident_cell() {
return ch;
@ -161,6 +163,7 @@ operator==(Triangulation_simplex_3<TriangulationDataStructure_3> s0,
typename Sim::Cell_handle neighbor;
switch (s0.dimension()) {
case -1: return s1.dimension() == -1;
case (0): // Vertex
return (s0.ch->vertex(s0.index(0)) == s1.ch->vertex(s1.index(0)));
case (1): // Edge
@ -180,7 +183,7 @@ operator==(Triangulation_simplex_3<TriangulationDataStructure_3> s0,
}
return false;
case (3):
return (&(*s0.ch) == &(*s1.ch));
return s0.ch.operator->() == s1.ch.operator->();
}
CGAL_error();
return false;

717
Triangulation_3/include/CGAL/Triangulation_3/internal/Triangulation_segment_traverser_3_impl.h
View File

@ -36,7 +36,7 @@ Triangulation_segment_cell_iterator_3( const Tr* tr, Vertex_handle s, Vertex_han
if( c->has_vertex( _tr->infinite_vertex(), inf ) )
c = c->neighbor(inf);
_cur = Simplex( c, Tr::VERTEX, c->index(s), -1 );
_cur = Simplex{ c, Tr::VERTEX, c->index(s), -1 };
jump_to_intersecting_cell();
}
@ -49,7 +49,7 @@ Triangulation_segment_cell_iterator_3( const Tr* tr, Vertex_handle s, const Poin
CGAL_precondition( s->point() != t );
CGAL_precondition( _tr->dimension() >= 2 );
CGAL_precondition( _tr->dimension() == 3 ||
orientation( *_tr->finite_facets_begin(), t ) == COPLANAR );
orientation( *_tr->finite_facets_begin(), t ) == COPLANAR );
_source = s->point();
_target = t;
@ -62,7 +62,7 @@ Triangulation_segment_cell_iterator_3( const Tr* tr, Vertex_handle s, const Poin
if( c->has_vertex( _tr->infinite_vertex(), inf ) )
c = c->neighbor(inf);
_cur = Simplex( c, Tr::VERTEX, c->index(s), -1 );
_cur = Simplex{ c, Tr::VERTEX, c->index(s), -1 };
jump_to_intersecting_cell();
}
@ -75,7 +75,7 @@ Triangulation_segment_cell_iterator_3( const Tr* tr, const Point& s, Vertex_hand
CGAL_precondition( s != t->point() );
CGAL_precondition( _tr->dimension() >= 2 );
CGAL_precondition( _tr->dimension() == 3 ||
orientation( *_tr->finite_facets_begin(), s ) == COPLANAR );
orientation( *_tr->finite_facets_begin(), s ) == COPLANAR );
_source = s;
_target = t->point();
@ -96,7 +96,7 @@ Triangulation_segment_cell_iterator_3( const Tr* tr, const Point& s, const Point
CGAL_precondition( s != t );
CGAL_precondition( _tr->dimension() >= 2 );
CGAL_precondition( _tr->dimension() == 3 ||
coplanar( *_tr->finite_facets_begin(), _target ) );
coplanar( *_tr->finite_facets_begin(), _target ) );
_source = s;
_target = t;
@ -124,7 +124,7 @@ template < class Tr, class Inc >
Triangulation_segment_cell_iterator_3<Tr,Inc>
Triangulation_segment_cell_iterator_3<Tr,Inc>::end() const {
SCI sci(_tr);
std::get<0>(sci._cur) = Cell_handle();
sci._cur.cell = Cell_handle();
return sci;
}
@ -239,8 +239,9 @@ walk_to_next() {
int ti;
if( cell()->has_vertex( _t_vertex, ti ) ) {
// The target is inside the cell.
_prev = Simplex( cell(), Tr::VERTEX, ti, -1 );
_prev = Simplex{ cell(), Tr::VERTEX, ti, -1 };
cell() = Cell_handle();
lt() = Locate_type::VERTEX;
return;
}
@ -302,24 +303,24 @@ bool Triangulation_segment_cell_iterator_3<Tr, Inc>::
have_same_entry(const Simplex& s1, const Simplex& s2) const
{
//type
if (std::get<1>(s1) != std::get<1>(s2))
if (s1.lt != s2.lt)
return false;
switch (std::get<1>(s1))
switch (s1.lt)
{
case Locate_type::VERTEX:
return std::get<0>(s1)->vertex(std::get<2>(s1)) == std::get<0>(s2)->vertex(std::get<2>(s2));
return s1.cell->vertex(s1.li) == s2.cell->vertex(s2.li);
case Locate_type::EDGE:
{
Vertex_handle v1a = std::get<0>(s1)->vertex(std::get<2>(s1));
Vertex_handle v1b = std::get<0>(s1)->vertex(std::get<3>(s1));
Vertex_handle v2a = std::get<0>(s2)->vertex(std::get<2>(s2));
Vertex_handle v2b = std::get<0>(s2)->vertex(std::get<3>(s2));
Vertex_handle v1a = s1.cell->vertex(s1.li);
Vertex_handle v1b = s1.cell->vertex(s1.lj);
Vertex_handle v2a = s2.cell->vertex(s2.li);
Vertex_handle v2b = s2.cell->vertex(s2.lj);
return (v1a == v2a && v1b == v2b)
|| (v1a == v2b && v1b == v2a);
}
case Locate_type::FACET:
return triangulation()->are_equal(Facet(std::get<0>(s1), std::get<2>(s1)),
Facet(std::get<0>(s2), std::get<2>(s2)));
return triangulation()->are_equal(Facet(s1.cell, s1.li),
Facet(s2.cell, s2.li));
default:
CGAL_assertion(false);
};
@ -332,305 +333,296 @@ std::pair<typename Triangulation_segment_cell_iterator_3<Tr, Inc>::Simplex,
Triangulation_segment_cell_iterator_3<Tr,Inc>::walk_to_next_3(const Simplex& prev,
const Simplex& cur) const
{
std::array<const Point*, 4> vert
= {&(std::get<0>(cur)->vertex(0)->point()),
&(std::get<0>(cur)->vertex(1)->point()),
&(std::get<0>(cur)->vertex(2)->point()),
&(std::get<0>(cur)->vertex(3)->point()) };
const auto cur_cell = cur.cell;
std::array<const Point*, 4> vert = {&(cur_cell->vertex(0)->point()), &(cur_cell->vertex(1)->point()),
&(cur_cell->vertex(2)->point()), &(cur_cell->vertex(3)->point())};
int inside=0,outside=0,regular_case=0,degenerate=0;
Cell_handle nnext;
int inside = 0, outside = 0, regular_case = 0, degenerate = 0;
if (std::get<1>(cur) == Tr::FACET) {
regular_case = 1;
int i = std::get<2>(cur);
int j0 = Tr::vertex_triple_index(i, 0);
int j1 = Tr::vertex_triple_index(i, 1);
int j2 = Tr::vertex_triple_index(i, 2);
Orientation o0 = _tr->orientation(_source, *vert[i], *vert[j0], _target);
if (o0 == POSITIVE) {
Orientation o1 = _tr->orientation(_source, *vert[i], *vert[j1], _target);
if (o1 != POSITIVE) {
if (_tr->orientation(*vert[i], *vert[j0], *vert[j1], _target) == POSITIVE) {
nnext = std::get<0>(cur)->neighbor(j2);
outside = j2;
if (o1 == ZERO) degenerate = 1; //EDGE i j1
}
else
inside = 1;
}
else {
if (_tr->orientation(*vert[i], *vert[j1], *vert[j2], _target) == POSITIVE) {
nnext = std::get<0>(cur)->neighbor(j0);
outside = j0;
}
else
inside = 2;
}
}
else if (o0 == ZERO) {
Orientation o1 = _tr->orientation(_source, *vert[i], *vert[j1], _target);
if (o1 == NEGATIVE) {
if (_tr->orientation(*vert[i], *vert[j0], *vert[j1], _target) == POSITIVE) {
nnext = std::get<0>(cur)->neighbor(j2); //EDGE i j0
degenerate = 2;
outside = 44;
}
else
inside = 3;
}
else if (o1 == ZERO) {
if (_tr->orientation(*vert[i], *vert[j0], *vert[j2], _target) == POSITIVE)
inside = 55;
else
{
nnext = std::get<0>(cur)->neighbor(j2); //VERTEX i
degenerate = 3;
outside = 5;
}
}
else {
if (_tr->orientation(*vert[i], *vert[j1], *vert[j2], _target) == POSITIVE) {
nnext = std::get<0>(cur)->neighbor(j0);
outside = j0;
}
else
inside = 4;
}
}
else {
Orientation o2 = _tr->orientation(_source, *vert[i], *vert[j2], _target);
if (o2 != NEGATIVE) {
if (_tr->orientation(*vert[i], *vert[j2], *vert[j0], _target) == POSITIVE) {
nnext = std::get<0>(cur)->neighbor(j1);
outside = j1;
if (o2 == ZERO) degenerate = 4; // EDGE i j2
}
else
inside = 5;
}
else {
if (_tr->orientation(*vert[i], *vert[j1], *vert[j2], _target) == POSITIVE) {
nnext = std::get<0>(cur)->neighbor(j0);
outside = j0;
}
else
inside = 6;
}
}
if ((!degenerate) && (!inside))
{
Simplex prev_after_walk(std::get<0>(cur), Tr::FACET, outside, -1);
Simplex cur_after_walk( nnext, Tr::FACET, nnext->index(std::get<0>(cur)), -1);
return std::make_pair(prev_after_walk, cur_after_walk);
}
if ((!degenerate) && inside)
{
Simplex prev_after_walk(std::get<0>(cur), Tr::CELL, -1, -1);
Simplex cur_after_walk(Cell_handle(), Tr::OUTSIDE_AFFINE_HULL, -1, -1);
return std::make_pair(prev_after_walk, cur_after_walk);
}
}
// We check in which direction the target lies
// by comparing its position relative to the planes through the
// source and the edges of the cell.
Orientation o[6];
Orientation op[4];
int pos = 0;
// We keep track of which orientations are calculated.
bool calc[6] = { false, false, false, false, false, false };
if( std::get<1>(cur) == Tr::VERTEX ) {
// The three planes through the vertex are set to coplanar.
for( int j = 0; j < 4; ++j ) {
if( std::get<2>(cur) != j ) {
int ij = edgeIndex( std::get<2>(cur), j );
o[ij] = COPLANAR;
calc[ij] = true;
}
}
}
else if( std::get<1>(cur) == Tr::EDGE ) {
// The plane through the edge is set to coplanar.
int ij = edgeIndex( std::get<2>(cur), std::get<3>(cur) );
o[ij] = COPLANAR;
calc[ij] = true;
}
// For the remembering stochastic walk, we start trying with a random facet.
int li = 0;
CGAL_assertion_code( bool incell = true; )
for( int k = 0; k < 4; ++k, ++li )
{
// Skip the previous cell.
Cell_handle next = std::get<0>(cur)->neighbor(li);
if( next == std::get<0>(prev) )
{
op[li] = POSITIVE;
pos += li;
continue;
}
const Point* backup = vert[li];
vert[li] = &_target;
// Check if the target is on the opposite side of the supporting plane.
op[li] = _tr->orientation( *vert[0], *vert[1], *vert[2], *vert[3] );
if( op[li] == POSITIVE )
pos += li;
if( op[li] != NEGATIVE ) {
vert[li] = backup;
continue;
}
CGAL_assertion_code( incell = false; )
// Check if the target is inside the 3-wedge with
// the source as apex and the facet as an intersection.
int lj = 0;
int Or = 0;
for( int l = 0; l < 4; ++l, ++lj ) {
if( li == lj )
continue;
// We check the orientation of the target compared to the plane
// Through the source and the edge opposite of ij.
int oij = 5 - edgeIndex( li, lj );
if( !calc[oij] ) {
const Point* backup2 = vert[lj];
vert[lj] = &_source;
o[oij] = _tr->orientation( *vert[0], *vert[1], *vert[2], *vert[3] );
vert[lj] = backup2;
calc[oij] = true;
}
if( o[oij] == POSITIVE ) {
// The target is not inside the pyramid.
// Invert the planes.
// This can be safely done because either
// they were not calculated yet,
// or they will no longer be used.
for( int j = 0; j < 4; ++j ) {
if( li == j ) continue;
int oij = 5 - edgeIndex( li, j );
o[oij] = -o[oij];
}
Or = 0;
break;
}
else
Or -= o[oij];
}
if( Or == 0 ) {
// Either the target is not inside the pyramid,
// or the pyramid is degenerate.
vert[li] = backup;
continue;
}
// The target is inside the pyramid.
Simplex prev_after_walk;
Simplex cur_after_walk;
std::get<0>(prev_after_walk) = std::get<0>(cur);
std::get<0>(cur_after_walk) = next;
switch( Or ) {
case 3:
std::get<1>(prev_after_walk) = Tr::FACET;
std::get<2>(prev_after_walk) = li;
std::get<1>(cur_after_walk) = Tr::FACET;
std::get<2>(cur_after_walk) = std::get<0>(cur_after_walk)->index(std::get<0>(prev_after_walk));
if(regular_case)
{
CGAL_assertion( std::get<0>(cur_after_walk)==nnext );
CGAL_assertion( li==outside );
CGAL_assertion( ! inside );
}
return std::make_pair(prev_after_walk, cur_after_walk);
case 2:
if(regular_case)
CGAL_assertion(degenerate );
std::get<1>(prev_after_walk) = Tr::EDGE;
std::get<1>(cur_after_walk) = Tr::EDGE;
for( int j = 0; j < 4; ++j ) {
if( li != j && o[ 5 - edgeIndex(li, j) ] == COPLANAR) {
Edge opp = opposite_edge( std::get<0>(prev), li, j );
std::get<2>(prev_after_walk) = opp.second;
std::get<3>(prev_after_walk) = opp.third;
std::get<2>(cur_after_walk)
= std::get<0>(cur_after_walk)->index(
std::get<0>(prev_after_walk)->vertex( std::get<2>(prev_after_walk) ) );
std::get<3>(cur_after_walk)
= std::get<0>(cur_after_walk)->index(
std::get<0>(prev_after_walk)->vertex( std::get<3>(prev_after_walk) ) );
return std::make_pair(prev_after_walk, cur_after_walk);
}
}
CGAL_assertion( false );
return std::make_pair(prev, cur);
case 1:
if(regular_case)
CGAL_assertion(degenerate );
std::get<1>(prev_after_walk) = Tr::VERTEX;
std::get<1>(cur_after_walk) = Tr::VERTEX;
for( int j = 0; j < 4; ++j ) {
if( li != j && o[ 5 - edgeIndex(li, j) ] == NEGATIVE ) {
std::get<2>(prev_after_walk) = j;
std::get<2>(cur_after_walk)
= std::get<0>(cur_after_walk)->index(
std::get<0>(prev_after_walk)->vertex(j) );
return std::make_pair(prev_after_walk, cur_after_walk);
}
}
CGAL_assertion( false );
return std::make_pair(prev, cur);
default:
CGAL_assertion( false );
return std::make_pair(prev, cur);
}
}
// The target lies inside this cell.
if(cur.lt == Tr::FACET && prev.cell != Cell_handle()) {
// [source, target] entered the cell `cur` via a facet.
// Note that, if prev.cell == Cell_handle(), that means `source` is *on*
// the facet, and the block of this `if` cannot be applied.
Simplex prev_after_walk;
CGAL_assertion( incell );
switch( op[0] + op[1] + op[2] + op[3] ) {
case 4:
CGAL_assertion( pos == 6 );
prev_after_walk = Simplex( std::get<0>(cur), Tr::CELL, -1, -1 );
CGAL_assertion( (! regular_case) || inside );
break;
Simplex cur_after_walk;
case 3:
prev_after_walk = Simplex( std::get<0>(cur), Tr::FACET, 6-pos, -1 );
break;
case 2:
if( pos < 3 )
prev_after_walk = Simplex( std::get<0>(cur), Tr::EDGE, 0, pos+1 );
else if( pos < 5 )
prev_after_walk = Simplex( std::get<0>(cur), Tr::EDGE, 1, pos-1 );
auto case_target_is_inside_cur_cell = [&](int case_nb) {
inside = case_nb;
prev_after_walk = {cur_cell, Tr::CELL, -1, -1};
cur_after_walk = {{}, Tr::CELL, -1, -1};
};
auto case_segment_exits_cur_cell_by = [&](int facet_nb, Cell_handle nnext = {}) {
if(nnext == Cell_handle{}) {
nnext = cur_cell->neighbor(facet_nb);
}
outside = facet_nb;
prev_after_walk = {cur_cell, Tr::FACET, facet_nb, -1};
cur_after_walk = {nnext, Tr::FACET, nnext->index(cur_cell), -1};
};
regular_case = 1;
const int i = cur.li;
const int j0 = Tr::vertex_triple_index(i, 0);
const int j1 = Tr::vertex_triple_index(i, 1);
const int j2 = Tr::vertex_triple_index(i, 2);
Orientation o0 = _tr->orientation(_source, *vert[i], *vert[j0], _target);
if(o0 == POSITIVE) { // o0 > 0
Orientation o1 = _tr->orientation(_source, *vert[i], *vert[j1], _target);
if(o1 != POSITIVE) { // o1 <= 0
Orientation oi01 = _tr->orientation(*vert[i], *vert[j0], *vert[j1], _target);
if(oi01 == POSITIVE) {
case_segment_exits_cur_cell_by(j2);
if(o1 == ZERO)
degenerate = 1; // EDGE i j1
} else { // o0 > 0, o1 <= 0, oi01 <= 0
case_target_is_inside_cur_cell(1);
if(oi01 == ZERO) { // on FACET j2 (i, j0, j1)
degenerate = 1;
} // end oi01 == ZERO
}
} // end o1 <= 0
else
prev_after_walk = Simplex( std::get<0>(cur), Tr::EDGE, 2, 3 );
break;
case 1:
prev_after_walk = Simplex( std::get<0>(cur), Tr::VERTEX, pos, -1 );
break;
default:
prev_after_walk = Simplex( std::get<0>(cur), Tr::OUTSIDE_AFFINE_HULL, -1, -1 );
CGAL_assertion( false );
{ // o1 > 0
Orientation oi12 = _tr->orientation(*vert[i], *vert[j1], *vert[j2], _target);
if(oi12 == POSITIVE) {
case_segment_exits_cur_cell_by(j0);
} else { // o0 > 0, o1 > 0, oi12 <= 0
case_target_is_inside_cur_cell(2);
if(oi12 == ZERO) { // on FACET j0 (i, j1, j2)
degenerate = 1;
} // end oi12 == ZERO
}
}
} // end o0 > 0
else if(o0 == ZERO)
{
// target is on plane (source, vert[i], vert[j0])
Orientation o1 = _tr->orientation(_source, *vert[i], *vert[j1], _target);
if(o1 == NEGATIVE) {
Orientation oi12 = _tr->orientation(*vert[i], *vert[j0], *vert[j1], _target);
if(oi12 == POSITIVE) {
degenerate = 2;
case_segment_exits_cur_cell_by(44, cur_cell->neighbor(j2)); // EDGE i j0
} else {
case_target_is_inside_cur_cell(3);
if(oi12 == ZERO) { // target is *on* EDGE i j0
degenerate = 1;
}
}
} else if(o1 == ZERO) {
// o0 == o1 == 0 -> target is on line source-vert[i]
if(_tr->orientation(*vert[i], *vert[j0], *vert[j2], _target) == POSITIVE)
case_target_is_inside_cur_cell(55);
else {
degenerate = 3;
case_segment_exits_cur_cell_by(5, cur_cell->neighbor(j2)); // VERTEX i
}
} else { // o0 == 0, o1 > 0
Orientation oi12 = _tr->orientation(*vert[i], *vert[j1], *vert[j2], _target);
if(oi12 == POSITIVE) {
case_segment_exits_cur_cell_by(j0);
} else {
case_target_is_inside_cur_cell(4);
if(oi12 == ZERO) { // on FACET j0 (i, j1, j2)
degenerate = 1;
} // end oi12 == ZERO
}
}
} // end o0 == 0
else
{ // o0 < 0
Orientation o2 = _tr->orientation(_source, *vert[i], *vert[j2], _target);
if(o2 != NEGATIVE) {
// o2 >= 0
Orientation oi20 = _tr->orientation(*vert[i], *vert[j2], *vert[j0], _target);
if(oi20 == POSITIVE) {
case_segment_exits_cur_cell_by(j1);
if(o2 == ZERO)
degenerate = 4; // EDGE i j2
} else {
case_target_is_inside_cur_cell(5);
if(oi20 == ZERO) { // on FACET j1 (i, j2, j0)
degenerate = 1;
}
}
} else {
Orientation oi12 = _tr->orientation(*vert[i], *vert[j1], *vert[j2], _target);
if(oi12 == POSITIVE) {
case_segment_exits_cur_cell_by(j0);
} else {
case_target_is_inside_cur_cell(6);
if(oi12 == ZERO) { // on FACET j0 (i, j1, j2)
degenerate = 1;
}
}
}
}
Simplex cur_after_walk(Cell_handle(), Tr::OUTSIDE_AFFINE_HULL, -1, -1);
return std::make_pair(prev_after_walk, cur_after_walk);
if(!degenerate) {
return {prev_after_walk, cur_after_walk};
}
}
// We check in which direction the target lies
// by comparing its position relative to the planes through the
// source and the edges of the cell.
std::array<Orientation, 6> o;
std::array<Orientation, 4> op;
int pos = 0;
// We keep track of which orientations are calculated.
bool calc[6] = {false, false, false, false, false, false};
if(cur.lt == Tr::VERTEX) {
// The three planes through the vertex are set to coplanar.
for(int j = 0; j < 4; ++j) {
if(cur.li != j) {
int ij = edgeIndex(cur.li, j);
o[ij] = COPLANAR;
calc[ij] = true;
}
}
} else if(cur.lt == Tr::EDGE) {
// The plane through the edge is set to coplanar.
int ij = edgeIndex(cur.li, cur.lj);
o[ij] = COPLANAR;
calc[ij] = true;
}
// For the remembering stochastic walk, we start trying with a random facet.
CGAL_assertion_code(bool incell = true;)
for(int li = 0; li < 4; ++li)
{
// Skip the previous cell.
Cell_handle next = cur_cell->neighbor(li);
if(next == prev.cell) {
op[li] = POSITIVE;
pos += li;
continue;
}
const Point* const backup_vert_li = std::exchange(vert[li], &_target);
// Check if the target is on the opposite side of the supporting plane.
op[li] = _tr->orientation(*vert[0], *vert[1], *vert[2], *vert[3]);
if(op[li] == POSITIVE)
pos += li;
if(op[li] != NEGATIVE) {
vert[li] = backup_vert_li;
continue;
}
CGAL_assertion_code(incell = false;)
// Check if the target is inside the 3-wedge with
// the source as apex and the facet as an intersection.
int Or = 0;
for(int lj = 0; lj < 4; ++lj) {
if(li == lj)
continue;
// We check the orientation of the target compared to the plane
// Through the source and the edge opposite of ij.
const int oij = 5 - edgeIndex(li, lj);
if(!calc[oij]) {
const Point* const backup_vert_lj = std::exchange(vert[lj], &_source);
o[oij] = _tr->orientation(*vert[0], *vert[1], *vert[2], *vert[3]);
vert[lj] = backup_vert_lj;
calc[oij] = true;
}
if(o[oij] == POSITIVE) {
// The target is not inside the pyramid.
// Invert the planes.
for(int j = 0; j < 4; ++j) {
if(li == j)
continue;
int oij = 5 - edgeIndex(li, j);
if(calc[oij])
o[oij] = -o[oij];
}
Or = 0;
break;
} else
Or -= o[oij];
}
if(Or == 0) {
// Either the target is not inside the pyramid,
// or the pyramid is degenerate.
vert[li] = backup_vert_li;
continue;
}
// The target is inside the pyramid.
switch(Or) {
case 3: {
if(regular_case) {
CGAL_assertion(li == outside);
CGAL_assertion(!inside);
}
return {{cur_cell, Tr::FACET, li}, {next, Tr::FACET, next->index(cur_cell)}};
}
case 2: {
if(regular_case)
CGAL_assertion(degenerate);
for(int j = 0; j < 4; ++j) {
if(li != j && o[5 - edgeIndex(li, j)] == COPLANAR) {
Edge opp = opposite_edge(prev.cell, li, j);
return {
{cur_cell, Tr::EDGE, opp.second, opp.third},
{next, Tr::EDGE, next->index(cur_cell->vertex(opp.second)), next->index(cur_cell->vertex(opp.third))}};
}
}
CGAL_unreachable();
return std::make_pair(prev, cur);
}
case 1:
if(regular_case)
CGAL_assertion(degenerate);
for(int j = 0; j < 4; ++j) {
if(li != j && o[5 - edgeIndex(li, j)] == NEGATIVE) {
return {{cur_cell, Tr::VERTEX, j}, {next, Tr::VERTEX, next->index(cur_cell->vertex(j))}};
}
}
CGAL_unreachable();
return std::make_pair(prev, cur);
default:
CGAL_unreachable();
return std::make_pair(prev, cur);
}
CGAL_unreachable();
}
// The target lies inside this cell.
CGAL_assertion( incell );
return {
[&]() -> Simplex {
switch( op[0] + op[1] + op[2] + op[3] ) {
case 4:
CGAL_assertion( pos == 6 );
CGAL_assertion( (! regular_case) || inside );
return { cur_cell, Tr::CELL };
break;
case 3:
return { cur_cell, Tr::FACET, 6 - pos };
break;
case 2:
if( pos < 3 ) // first is 0
return { cur_cell, Tr::EDGE, 0, pos };
else if( pos < 5 ) { // could be (0, pos), or (1, pos-1)
if(op[0] == POSITIVE)
return { cur_cell, Tr::EDGE, 0, pos };
else
return { cur_cell, Tr::EDGE, 1, pos-1 };
}
else
return { cur_cell, Tr::EDGE, 2, 3 };
break;
case 1:
return { cur_cell, Tr::VERTEX, pos };
break;
default:
CGAL_unreachable();
}
}(),
{ Cell_handle() }
};
}
template < class Tr, class Inc >
@ -643,7 +635,9 @@ walk_to_next_3_inf( int inf )
Cell_handle fin = cell()->neighbor(inf);
if( fin == prev_cell() ) {
_prev = _cur;
prev_lt() = Tr::CELL;
cell() = Cell_handle();
lt() = Tr::CELL;
return;
}
@ -658,7 +652,7 @@ walk_to_next_3_inf( int inf )
if( _tr->orientation( *vert[0], *vert[1], *vert[2], *vert[3] ) == POSITIVE ) {
// The target lies in an infinite cell.
// Note that we do not traverse to other infinite cells.
_prev = Simplex( cell(), Tr::OUTSIDE_CONVEX_HULL, -1, -1 );
_prev = Simplex{ cell(), Tr::OUTSIDE_CONVEX_HULL, -1, -1 };
cell() = Cell_handle();
return;
}
@ -683,20 +677,20 @@ walk_to_next_3_inf( int inf )
continue;
}
Point* backup = vert[li];
Point* backup_vert_li = vert[li];
vert[li] = &(_target);
o[li] = _tr->orientation( *vert[0], *vert[1], *vert[2], *vert[3] );
if( o[li] != NEGATIVE ) {
vert[li] = backup;
vert[li] = backup_vert_li;
continue;
}
// The target lies behind the plane through the source and two finite vertices.
// Traverse to the incident infinite cell.
CGAL_assertion( _tr->is_infinite( next ) );
_prev = Simplex( cell(), Tr::FACET, li, -1 );
_cur = Simplex( next, Tr::FACET, next->index( prev_cell() ), -1 );
_prev = Simplex{ cell(), Tr::FACET, li, -1 };
_cur = Simplex{ next, Tr::FACET, next->index(prev_cell()), -1 };
return;
}
@ -725,7 +719,7 @@ walk_to_next_3_inf( int inf )
return;
}
}
CGAL_assertion( false );
CGAL_unreachable();
return;
case 1:
prev_lt() = Tr::VERTEX;
@ -737,10 +731,10 @@ walk_to_next_3_inf( int inf )
return;
}
}
CGAL_assertion( false );
CGAL_unreachable();
return;
default:
CGAL_assertion( false );
CGAL_unreachable();
return;
}
}
@ -802,7 +796,7 @@ walk_to_next_2()
return;
default:
// The current vertex is the target.
CGAL_assertion(false);
CGAL_unreachable();
return;
}
}
@ -810,27 +804,26 @@ walk_to_next_2()
// The target lies in this cell.
switch( ocw+occw+op ) {
case 3:
_prev = Simplex( cell(), Tr::FACET, 3, -1 );
break;
_prev = Simplex{ cell(), Tr::FACET, 3, -1 };
break;
case 2:
if( ocw == 0 )
_prev = Simplex( cell(), Tr::EDGE, _tr->ccw(li()), -1 );
_prev = Simplex{ cell(), Tr::EDGE, _tr->ccw(li()), -1 };
else if( occw == 0 )
_prev = Simplex( cell(), Tr::EDGE, _tr->cw(li()), -1 );
_prev = Simplex{ cell(), Tr::EDGE, _tr->cw(li()), -1 };
else
_prev = Simplex( cell(), Tr::EDGE, li(), -1 );
_prev = Simplex{ cell(), Tr::EDGE, li(), -1 };
break;
case 1:
if( ocw == 1 )
_prev = Simplex( cell(), Tr::VERTEX, _tr->ccw(li()), -1 );
_prev = Simplex{ cell(), Tr::VERTEX, _tr->ccw(li()), -1 };
else if( occw == 1 )
_prev = Simplex( cell(), Tr::VERTEX, _tr->cw(li()), -1 );
_prev = Simplex{ cell(), Tr::VERTEX, _tr->cw(li()), -1 };
else
_prev = Simplex( cell(), Tr::VERTEX, li(), -1 );
_prev = Simplex{ cell(), Tr::VERTEX, li(), -1 };
break;
case 0:
CGAL_assertion(false);
_prev = Simplex( cell(), Tr::OUTSIDE_AFFINE_HULL, -1, -1 );
CGAL_unreachable();
break;
}
cell() = Cell_handle();
@ -927,19 +920,19 @@ walk_to_next_2()
// The target lies in this cell.
if( op == POSITIVE )
_prev = Simplex( cell(), Tr::FACET, 3, -1 );
_prev = Simplex{ cell(), Tr::FACET, 3, -1 };
else {
CGAL_assertion( op == ZERO );
switch( o ) {
case POSITIVE:
_prev = Simplex( cell(), Tr::EDGE, li(), lk );
break;
_prev = Simplex{ cell(), Tr::EDGE, li(), lk };
break;
case NEGATIVE:
_prev = Simplex( cell(), Tr::EDGE, lj(), lk );
break;
_prev = Simplex{ cell(), Tr::EDGE, lj(), lk };
break;
case ZERO:
_prev = Simplex( cell(), Tr::VERTEX, lk, -1 );
break;
_prev = Simplex{ cell(), Tr::VERTEX, lk, -1 };
break;
}
}
cell() = Cell_handle();
@ -969,7 +962,7 @@ walk_to_next_2()
if( o[_tr->ccw(li)] == NEGATIVE )
continue;
else if( op == COLLINEAR && o[_tr->ccw(li)] == COLLINEAR ) {
_prev = Simplex( cell(), Tr::VERTEX, _tr->ccw(li), -1 );
_prev = Simplex{ cell(), Tr::VERTEX, _tr->ccw(li), -1 };
cell() = Cell_handle();
return;
}
@ -981,7 +974,7 @@ walk_to_next_2()
if( o[_tr->cw(li)] == POSITIVE )
continue;
else if( op == COLLINEAR && o[_tr->cw(li)] == COLLINEAR ) {
_prev = Simplex( cell(), Tr::VERTEX, _tr->cw(li), -1 );
_prev = Simplex{ cell(), Tr::VERTEX, _tr->cw(li), -1 };
cell() = Cell_handle();
return;
}
@ -1006,18 +999,18 @@ walk_to_next_2()
this->li() = cell()->index( prev_cell()->vertex( prev_li() ) );
return;
default:
CGAL_assertion( false );
CGAL_unreachable();
return;
}
}
// The target lies in this cell.
_prev = Simplex( cell(), Tr::FACET, 3, -1 );
_prev = Simplex{ cell(), Tr::FACET, 3, -1 };
cell() = Cell_handle();
return;
}
default:
CGAL_assertion( false );
CGAL_unreachable();
}
}
@ -1046,8 +1039,8 @@ walk_to_next_2_inf( int inf )
_target );
if( occw == NEGATIVE ) {
Cell_handle tmp = cell()->neighbor(_tr->cw(inf));
_prev = Simplex( cell(), Tr::EDGE, _tr->ccw(inf), inf );
_cur = Simplex( tmp, Tr::EDGE, tmp->index( prev_cell()->vertex( prev_li() ) ), tmp->index( prev_cell()->vertex( prev_lj() ) ) );
_prev = Simplex{ cell(), Tr::EDGE, _tr->ccw(inf), inf };
_cur = Simplex{ tmp, Tr::EDGE, tmp->index(prev_cell()->vertex(prev_li())), tmp->index(prev_cell()->vertex(prev_lj())) };
return;
}
Orientation ocw = coplanar_orientation( _source,
@ -1056,8 +1049,8 @@ walk_to_next_2_inf( int inf )
_target );
if( ocw == NEGATIVE ) {
Cell_handle tmp = cell()->neighbor(_tr->ccw(inf));
_prev = Simplex( cell(), Tr::EDGE, _tr->cw(inf), inf );
_cur = Simplex( tmp, Tr::EDGE, tmp->index( prev_cell()->vertex( prev_li() ) ), tmp->index( prev_cell()->vertex( prev_lj() ) ) );
_prev = Simplex{ cell(), Tr::EDGE, _tr->cw(inf), inf };
_cur = Simplex{ tmp, Tr::EDGE, tmp->index(prev_cell()->vertex(prev_li())), tmp->index(prev_cell()->vertex(prev_lj())) };
return;
}
Orientation op = coplanar_orientation(
@ -1067,37 +1060,37 @@ walk_to_next_2_inf( int inf )
switch( op ) {
case NEGATIVE:
if( occw == COLLINEAR ) {
_prev = Simplex( cell(), Tr::VERTEX, _tr->ccw(inf), -1 );
_cur = Simplex( fin, Tr::VERTEX, fin->index( prev_cell()->vertex( prev_li() ) ), -1 );
return;
_prev = Simplex{ cell(), Tr::VERTEX, _tr->ccw(inf), -1 };
_cur = Simplex{ fin, Tr::VERTEX, fin->index(prev_cell()->vertex(prev_li())), -1 };
return;
}
if( ocw == COLLINEAR ) {
_prev = Simplex( cell(), Tr::VERTEX, _tr->cw(inf), -1 );
_cur = Simplex( fin, Tr::VERTEX, fin->index( prev_cell()->vertex( prev_li() ) ), -1 );
return;
_prev = Simplex{ cell(), Tr::VERTEX, _tr->cw(inf), -1 };
_cur = Simplex{ fin, Tr::VERTEX, fin->index(prev_cell()->vertex(prev_li())), -1 };
return;
}
_prev = Simplex( cell(), Tr::EDGE, _tr->ccw(inf), _tr->cw(inf) );
_cur = Simplex( fin, Tr::EDGE, fin->index( prev_cell()->vertex( prev_li() ) ), fin->index( prev_cell()->vertex( prev_lj() ) ) );
_prev = Simplex{ cell(), Tr::EDGE, _tr->ccw(inf), _tr->cw(inf) };
_cur = Simplex{ fin, Tr::EDGE, fin->index(prev_cell()->vertex(prev_li())), fin->index(prev_cell()->vertex(prev_lj())) };
return;
case COLLINEAR:
if( occw == COLLINEAR ) {
_prev = Simplex( cell(), Tr::VERTEX, _tr->ccw(inf), -1 );
cell() = Cell_handle();
return;
_prev = Simplex{ cell(), Tr::VERTEX, _tr->ccw(inf), -1 };
cell() = Cell_handle();
return;
}
if( ocw == COLLINEAR ) {
_prev = Simplex( cell(), Tr::VERTEX, _tr->cw(inf), -1 );
cell() = Cell_handle();
return;
_prev = Simplex{ cell(), Tr::VERTEX, _tr->cw(inf), -1 };
cell() = Cell_handle();
return;
}
_prev = Simplex( cell(), Tr::EDGE, _tr->ccw(inf), _tr->cw(inf) );
_prev = Simplex{ cell(), Tr::EDGE, _tr->ccw(inf), _tr->cw(inf) };
cell() = Cell_handle();
return;
case POSITIVE:
// The tarstd::std::get lies in this infinite cell.
_prev = Simplex( cell(), Tr::OUTSIDE_CONVEX_HULL, -1, -1 );
cell() = Cell_handle();
return;
_prev = Simplex{ cell(), Tr::OUTSIDE_CONVEX_HULL, -1, -1 };
cell() = Cell_handle();
return;
}
}
@ -1138,7 +1131,7 @@ Triangulation_segment_cell_iterator_3<Tr, Inc>::opposite_edge(
case 5: return Edge(c, 2, 3);
}
CGAL_assertion(false);
CGAL_unreachable();
return Edge();
}

558
Triangulation_3/include/CGAL/Triangulation_segment_traverser_3.h
View File

@ -17,7 +17,6 @@
#include <iostream>
#include <utility>
#include <tuple>
#include <CGAL/assertions.h>
#include <CGAL/Triangulation_utils_3.h>
@ -27,6 +26,7 @@
#include <CGAL/Triangulation_vertex_base_3.h>
#include <CGAL/Triangulation_simplex_3.h>
#include <boost/optional.hpp>
// If defined, type casting is done statically,
// reducing type-safety overhead.
@ -113,7 +113,13 @@ public:
typedef typename Tr::Locate_type Locate_type; //< defines the simplex type returned from location.
typedef std::tuple<Cell_handle,Locate_type,int,int> Simplex; //< defines the simplex type.
struct Simplex //< defines the simplex type
{
Cell_handle cell = {};
Locate_type lt = Locate_type::OUTSIDE_AFFINE_HULL;
int li = -1;
int lj = -1;
};
typedef Cell value_type; //< defines the value type the iterator refers to.
typedef Cell& reference; //< defines the reference type of the iterator.
@ -127,6 +133,55 @@ public:
template < class Tr2, class Inc2 >
struct Rebind { typedef Triangulation_segment_cell_iterator_3<Tr2,Inc2> Other; };
#if CGAL_DEBUG_TRIANGULATION_SEGMENT_TRAVERSER_3
static auto display_vert(Vertex_handle v)
{
std::stringstream os;
os.precision(17);
if(v->time_stamp() == 0) {
os << "inf";
} else {
os << '#' << v->time_stamp() << "=(" << v->point() << ")";
}
return os.str();
};
static auto display_lt(Locate_type lt) {
std::stringstream os;
switch(lt) {
case Locate_type::VERTEX: os << " VERTEX"; break;
case Locate_type::EDGE: os << " EDGE"; break;
case Locate_type::FACET: os << " FACET"; break;
case Locate_type::CELL: os << " CELL"; break;
case Locate_type::OUTSIDE_CONVEX_HULL: os << " OUTSIDE_CONVEX_HULL"; break;
case Locate_type::OUTSIDE_AFFINE_HULL: os << " OUTSIDE_AFFINE_HULL"; break;
}
return os.str();
}
static auto debug_simplex(Simplex s) {
std::stringstream os;
os.precision(17);
const auto [c, lt, i, j] = s;
if(c == Cell_handle{}) {
os << "end()";
} else {
os << display_vert(c->vertex(0)) << " - " << display_vert(c->vertex(1)) << " - "
<< display_vert(c->vertex(2)) << " - " << display_vert(c->vertex(3));
os << display_lt(lt) << " " << i << " " << j;
}
return os.str();
}
auto debug_iterator() const
{
std::stringstream os;
os.precision(17);
os << " prev: " << debug_simplex(_prev) << "\n cur: " << debug_simplex(_cur);
return os.str();
}
#endif // CGAL_DEBUG_TRIANGULATION_SEGMENT_TRAVERSER_3
private:
typedef Segment_cell_iterator SCI;
@ -228,15 +283,17 @@ public:
*/
const Point& target() const { return _target; }
Vertex_handle target_vertex() const { return _t_vertex; }
// gives a handle to the current cell.
/* By invariance, this cell is intersected by the segment
* between `source()` and `target()`.
* \return a handle to the current cell.
* \sa `cell()`.
*/
Cell_handle handle()
Cell_handle handle() const
{
return std::get<0>(_cur);
return _cur.cell;
}
// gives the previous cell.
@ -257,7 +314,7 @@ public:
*/
Cell* operator->()
{
return &*std::get<0>(_cur);
return &*(_cur.cell);
}
// provides an indirection operator.
@ -265,7 +322,7 @@ public:
*/
Cell& operator*()
{
return *std::get<0>(_cur);
return *(_cur.cell);
}
// provides a conversion operator.
@ -273,7 +330,7 @@ public:
*/
operator Cell_handle() const
{
return std::get<0>(_cur);
return _cur.cell;
}
// provides a conversion operator.
@ -299,6 +356,10 @@ public:
{
lt = this->lt(); li = this->li(); lj = this->lj();
}
std::tuple<Locate_type, int, int> entry() const
{
return { lt(), li(), lj() };
}
// gives the simplex through which the previous cell was exited.
/* \pre the current cell is not the initial cell.
*/
@ -306,6 +367,10 @@ public:
{
lt = prev_lt(); li = prev_li(); lj = prev_lj();
}
std::tuple<Locate_type, int, int> exit() const
{
return { prev_lt(), prev_li(), prev_lj() };
}
// gives the past-the-end iterator associated with this iterator.
SCI end() const;
@ -364,7 +429,7 @@ public:
*/
bool operator==( const Cell_handle& ch ) const
{
return ch == std::get<0>(_cur);
return ch == _cur.cell;
}
// compares the current cell with `ch`.
@ -375,7 +440,7 @@ public:
*/
bool operator!=( const Cell_handle& ch ) const
{
return ch != std::get<0>(_cur);
return ch != _cur.cell;
}
// \}
@ -439,32 +504,33 @@ private:
Edge opposite_edge(Cell_handle c, int li, int lj) const;
Edge opposite_edge(const Edge& e) const;
protected:
// ref-accessors to the simplex, for use in internal code
// access _cur
Cell_handle& cell() { return std::get<0>(_cur); }
Cell_handle const& cell() const { return std::get<0>(_cur); }
Cell_handle& cell() { return _cur.cell; }
Cell_handle const& cell() const { return _cur.cell; }
Locate_type& lt() { return std::get<1>(_cur); }
Locate_type const& lt() const { return std::get<1>(_cur); }
Locate_type& lt() { return _cur.lt; }
Locate_type const& lt() const { return _cur.lt; }
int& li() { return std::get<2>(_cur); }
int const& li() const { return std::get<2>(_cur); }
int& li() { return _cur.li; }
int const& li() const { return _cur.li; }
int& lj() { return std::get<3>(_cur); }
int const& lj() const { return std::get<3>(_cur); }
int& lj() { return _cur.lj; }
int const& lj() const { return _cur.lj; }
// access _prev
Cell_handle& prev_cell() { return std::get<0>(_prev); }
Cell_handle const& prev_cell() const { return std::get<0>(_prev); }
Cell_handle& prev_cell() { return _prev.cell; }
Cell_handle const& prev_cell() const { return _prev.cell; }
Locate_type& prev_lt() { return std::get<1>(_prev); }
Locate_type const& prev_lt() const { return std::get<1>(_prev); }
Locate_type& prev_lt() { return _prev.lt; }
Locate_type const& prev_lt() const { return _prev.lt; }
int& prev_li() { return std::get<2>(_prev); }
int const& prev_li() const { return std::get<2>(_prev); }
int& prev_li() { return _prev.li; }
int const& prev_li() const { return _prev.li; }
int& prev_lj() { return std::get<3>(_prev); }
int const& prev_lj() const { return std::get<3>(_prev); }
int& prev_lj() { return _prev.lj; }
int const& prev_lj() const { return _prev.lj; }
}; // class Triangulation_segment_cell_iterator_3
@ -572,7 +638,7 @@ public:
const Point& source() const { return _cell_iterator.source(); }
const Point& target() const { return _cell_iterator.target(); }
const Tr* triangulation() const { return _cell_iterator.triangulation(); }
const Tr& triangulation() const { return *_cell_iterator.triangulation(); }
private:
Triangulation_segment_simplex_iterator_3
@ -584,20 +650,23 @@ private:
private:
void set_curr_simplex_to_entry()
{
#if CGAL_DEBUG_TRIANGULATION_SEGMENT_TRAVERSER_3
std::cerr << "cell iterator is:\n" << _cell_iterator.debug_iterator() << std::endl;
#endif // #if CGAL_DEBUG_TRIANGULATION_SEGMENT_TRAVERSER_3
Locate_type lt;
int li, lj;
Cell_handle cell;
Cell_handle cell = Cell_handle(_cell_iterator);
//check what is the entry type of _cell_iterator
if (Cell_handle(_cell_iterator) == Cell_handle())
if (cell == Cell_handle())
{
//where did the segment std::get out from previous cell
//where did the segment get out from previous cell
cell = _cell_iterator.previous();
_cell_iterator.exit(lt, li, lj);
}
else
{
cell = Cell_handle(_cell_iterator);
_cell_iterator.entry(lt, li, lj);
}
@ -623,7 +692,7 @@ private:
_curr_simplex = cell;
break;
default:
CGAL_assertion(false);
CGAL_unreachable();
};
}
@ -637,282 +706,145 @@ public:
// provides the increment postfix operator.
Simplex_iterator& operator++()
{
auto increment_cell_iterator = [&]() {
++_cell_iterator;
#if CGAL_DEBUG_TRIANGULATION_SEGMENT_TRAVERSER_3
std::cerr << "increment cell iterator to:\n" << _cell_iterator.debug_iterator() << '\n';
#endif
};
CGAL_assertion(_curr_simplex.incident_cell() != Cell_handle());
switch(_curr_simplex.dimension())
{
case 3 :/*Cell_handle*/
{
Cell_handle ch = Cell_handle(_cell_iterator);
if (ch == Cell_handle())
{
if(!triangulation()->is_infinite(Cell_handle(_curr_simplex)))
set_curr_simplex_to_entry();
else
_curr_simplex = Simplex_3();
break;
}
else
{
if (!cell_iterator_is_ahead())
++_cell_iterator;
set_curr_simplex_to_entry();
}
break;
if(!cell_iterator_is_ahead()) {
increment_cell_iterator(); // cell_iterator needs to be ahead
}
case 2 :/*Facet*/
Cell_handle ch_next = Cell_handle(_cell_iterator);
Cell_handle ch_prev = _cell_iterator.previous();
Locate_type lt_prev;
int li_prev, lj_prev;
_cell_iterator.exit(lt_prev, li_prev, lj_prev);
if(_curr_simplex.dimension() == 3) {
set_curr_simplex_to_entry();
return *this;
}
if(lt_prev == Locate_type::CELL ||
lt_prev == Locate_type::OUTSIDE_CONVEX_HULL ||
lt_prev == Locate_type::OUTSIDE_AFFINE_HULL)
{
Cell_handle ch = Cell_handle(_cell_iterator);
if (!cell_iterator_is_ahead())
{
//cell_iterator is not ahead. get_facet() is part of cell_iterator
//we cannot be in any of the degenerate cases, only detected by
//taking cell_iterator one step forward
CGAL_assertion(cell_has_facet(Cell_handle(_cell_iterator), get_facet()));
++_cell_iterator;
if (Cell_handle(_cell_iterator) == Cell_handle())
{
_curr_simplex = _cell_iterator.previous();
break;
}
}
else
ch = _cell_iterator.previous();
CGAL_assertion(ch_next == Cell_handle());
_curr_simplex = ch_prev;
return *this;
}
Cell_handle chnext = Cell_handle(_cell_iterator);
Locate_type ltnext;
int linext, ljnext;
_cell_iterator.entry(ltnext, linext, ljnext);
switch (ltnext)//entry simplex in next cell
{
case Locate_type::VERTEX:
{
if (_cell_iterator == _cell_iterator.end())
{
_curr_simplex = Simplex_3();
break;
}
//if the entry vertex is a vertex of current facet
int i;
if (triangulation()->has_vertex(get_facet(), chnext->vertex(linext), i))
set_curr_simplex_to_entry();
switch(_curr_simplex.dimension()) {
case 2: { /*Facet*/
CGAL_assertion((ch_next == Cell_handle()) == (_cell_iterator == _cell_iterator.end()));
switch(lt_prev) {
case Locate_type::VERTEX: { // facet-cell?-vertex-outside
Vertex_handle v_prev{ch_prev->vertex(li_prev)};
if(facet_has_vertex(get_facet(), v_prev))
_curr_simplex = v_prev;
else
_curr_simplex = ch;
break;
}
case Locate_type::EDGE:
if (facet_has_edge(get_facet(), Edge(chnext, linext, ljnext)))
set_curr_simplex_to_entry();
_curr_simplex = ch_prev;
} break;
case Locate_type::EDGE: { // facet-cell?-edge-outside
Edge edge_prev{ch_prev, li_prev, lj_prev};
if(facet_has_edge(get_facet(), edge_prev))
_curr_simplex = edge_prev;
else
_curr_simplex = ch;
break;
case Locate_type::FACET:
_curr_simplex = ch;
break;
case Locate_type::OUTSIDE_AFFINE_HULL:
_curr_simplex = Simplex_3();
break;
_curr_simplex = ch_prev;
} break;
case Locate_type::FACET: { // facet-cell-facet-outside
Facet f_prev{ch_prev, li_prev};
if(is_same_facet(f_prev, get_facet())) {
if(ch_next == Cell_handle())
_curr_simplex = Simplex_3();
else
_curr_simplex = ch_next;
} else
_curr_simplex = ch_prev;
} break;
default:
CGAL_assertion(false);
};
break;
}
case 1:/*Edge*/
{
Cell_handle ch = Cell_handle(_cell_iterator);
if (ch == _cell_iterator.previous())
{
_curr_simplex = Simplex_3();
break;
CGAL_unreachable();
}
Locate_type lt;
int li, lj;
_cell_iterator.entry(lt, li, lj);
if (!cell_iterator_is_ahead())
{
++_cell_iterator;//cell_iterator needs to be ahead to detect degeneracies
if (Cell_handle(_cell_iterator) == Cell_handle())
{
_curr_simplex = _cell_iterator.previous();
break;
} break;
case 1: {/*Edge*/
switch(lt_prev) {
case Locate_type::VERTEX: { //edge-vertex-outside
Vertex_handle v_prev{ch_prev->vertex(li_prev)};
if(edge_has_vertex(get_edge(), v_prev))
_curr_simplex = v_prev;
else
_curr_simplex = shared_facet(get_edge(), v_prev);
} break;
case Locate_type::EDGE: { //edge-outside or edge-cell-edge-outside
const Edge e_prev(ch_prev, li_prev, lj_prev);
if(is_same_edge(get_edge(), e_prev)) {
if(ch_next == Cell_handle()) {
_curr_simplex = Simplex_3();
} else {
_curr_simplex = ch_next;
}
} else {
auto facet_opt = shared_facet(get_edge(), e_prev);
if(static_cast<bool>(facet_opt)) {
_curr_simplex = *facet_opt;
}
else {
_curr_simplex = shared_cell(get_edge(), e_prev);
}
}
}
Cell_handle chnext = Cell_handle(_cell_iterator);
Locate_type ltnext;
int linext, ljnext;
_cell_iterator.entry(ltnext, linext, ljnext);
switch (ltnext)//entry simplex in next cell
{
case Locate_type::VERTEX:
if (edge_has_vertex(get_edge(), chnext->vertex(linext)))
_curr_simplex = chnext->vertex(linext);
} break;
case Locate_type::FACET: {
Facet f_prev{ch_prev, li_prev};
if(facet_has_edge(f_prev, get_edge()))
_curr_simplex = f_prev; //edge-facet-outside
else
_curr_simplex = shared_facet(get_edge(), chnext->vertex(linext));
break;
case Locate_type::EDGE:
{
CGAL_assertion(_cell_iterator == _cell_iterator.end()
|| triangulation()->is_infinite(chnext)
|| _curr_simplex != Simplex_3(Edge(chnext, linext, ljnext)));
if (_cell_iterator == _cell_iterator.end())
_curr_simplex = Simplex_3();
else if (triangulation()->is_infinite(chnext)
&& _curr_simplex == Simplex_3(Edge(chnext, linext, ljnext)))
_curr_simplex = chnext;
else
_curr_simplex = shared_facet(get_edge(), Edge(chnext, linext, ljnext));
break;
}
case Locate_type::FACET:
_curr_simplex = Cell_handle(_cell_iterator);//query goes through the cell
break;
case Locate_type::OUTSIDE_AFFINE_HULL:
{
Cell_handle chprev = _cell_iterator.previous();
Locate_type ltprev;
int liprev, ljprev;
_cell_iterator.exit(ltprev, liprev, ljprev);
if (ltprev == Locate_type::VERTEX) //edge-vertex-outside
_curr_simplex = chprev->vertex(liprev);
else
_curr_simplex = Simplex_3(); //edge-outside
break;
}
_curr_simplex = ch_prev; //query goes through the cell
} break;
default:
CGAL_assertion(false);//should not happen
};
break;
}
CGAL_unreachable();
}
} break;
case 0 :/*Vertex_handle*/
{
Cell_handle ch = Cell_handle(_cell_iterator);
if (ch == _cell_iterator.previous())
{
_curr_simplex = Simplex_3();
break;
}
if (!cell_iterator_is_ahead()) //_curr_simplex does contain v
{
++_cell_iterator;//cell_iterator needs to be ahead to detect degeneracies
}
else
ch = _cell_iterator.previous();
Cell_handle chnext = Cell_handle(_cell_iterator);
//_cell_iterator is one step forward _curr_simplex
CGAL_assertion(ch != chnext);
Locate_type ltnext;
int linext, ljnext;
_cell_iterator.entry(ltnext, linext, ljnext);
Cell_handle prev;
Locate_type ltprev;
int liprev, ljprev;
prev = _cell_iterator.previous();
_cell_iterator.exit(ltprev, liprev, ljprev);
switch (ltnext)
{
case Locate_type::VERTEX:
{
CGAL_assertion(_cell_iterator == _cell_iterator.end()
|| get_vertex() != chnext->vertex(linext)
|| triangulation()->is_infinite(chnext));
if (_cell_iterator == _cell_iterator.end())
{
if (prev == ch && ltprev == Locate_type::VERTEX)
{
CGAL_assertion(prev->vertex(liprev) == get_vertex());
_curr_simplex = ch;
}
else
{
if(ltprev == Locate_type::FACET)
_curr_simplex = Facet(prev, liprev);
else if(ltprev == Locate_type::EDGE)
_curr_simplex = Edge(prev, liprev, ljprev);
else
CGAL_assertion(false);
}
}
else
{
if (triangulation()->is_infinite(chnext) && get_vertex() == chnext->vertex(linext))
_curr_simplex = chnext;
else
{
Cell_handle ec;
int ei = -1, ej = -1;
if (!triangulation()->is_edge(get_vertex(), chnext->vertex(linext), ec, ei, ej))
CGAL_assertion(false);
_curr_simplex = Edge(ec, ei, ej);
}
}
break;
}
case Locate_type::EDGE:
{
//facet shared by get_vertex() and the edge
//none of ch and chnext is certainly shared by both endpoints
_curr_simplex = shared_facet(Edge(chnext, linext, ljnext), get_vertex());
break;
}
case Locate_type::OUTSIDE_AFFINE_HULL:
{
CGAL_assertion(_cell_iterator == _cell_iterator.end());
if (ltprev == Locate_type::VERTEX) //vertex-edge-vertex-outside
{
if(prev->vertex(liprev) != get_vertex())//avoid infinite loop edge-vertex-same edge-...
_curr_simplex = Edge(prev, liprev, prev->index(get_vertex()));
else
switch(lt_prev) {
case Locate_type::VERTEX: {
if(ch_prev->vertex(li_prev) != get_vertex()) // avoid infinite loop edge-vertex-same edge-...
_curr_simplex = Edge(ch_prev, li_prev, ch_prev->index(get_vertex()));
else {
if(ch_next == Cell_handle()) {
_curr_simplex = Simplex_3();
} else {
_curr_simplex = ch_next;
}
}
else if (ltprev == Locate_type::EDGE)//vertex-facet-edge-outside
_curr_simplex = Facet(prev, prev->index(get_vertex()));
else if (ltprev == Locate_type::FACET) //vertex-facet-outside
{
if(prev->vertex(liprev) != get_vertex()) //vertex-facet-outside
_curr_simplex = Facet(prev, liprev);
else //vertex-cell-facet-outside
_curr_simplex = prev;
}
} break;
case Locate_type::EDGE: {
const Edge e_prev(ch_prev, li_prev, lj_prev);
if(edge_has_vertex(e_prev, get_vertex()))
_curr_simplex = e_prev;
else
{
CGAL_assertion(ltprev == Locate_type::CELL);//vertex-cell-outside
_curr_simplex = prev;
}
break;
_curr_simplex = shared_facet(Edge(ch_prev, li_prev, lj_prev), get_vertex());
} break;
case Locate_type::FACET: {
if(ch_prev->vertex(li_prev) != get_vertex()) // vertex-facet-outside
_curr_simplex = Facet(ch_prev, li_prev);
else // vertex-cell-facet-outside
_curr_simplex = ch_prev;
} break;
default:
CGAL_unreachable();
}
default://FACET
if (chnext == Cell_handle())
_curr_simplex = Simplex_3();
else
_curr_simplex = shared_cell(Facet(chnext, linext), get_vertex());
break;
};
}
break;
} break;
default:
CGAL_assertion(false);
CGAL_unreachable();
};
return *this;
}
// provides the increment prefix operator.
Simplex_iterator operator++(int)
{
@ -1010,10 +942,10 @@ private:
case 3 ://cell
return ch != get_cell();
default:
CGAL_assertion(false);
CGAL_unreachable();
}
//should not be reached
CGAL_assertion(false);
CGAL_unreachable();
return false;
}
@ -1048,13 +980,29 @@ private:
return false;
}
bool facet_has_vertex(const Facet& f, const Vertex_handle v) const
{
return triangulation().tds().has_vertex(f, v);
}
bool edge_has_vertex(const Edge& e, const Vertex_handle v) const
{
return e.first->vertex(e.second) == v
|| e.first->vertex(e.third) == v;
}
Vertex_handle shared_vertex(const Edge& e1, const Edge& e2) const
bool is_same_edge(const Edge& e1, const Edge& e2) const
{
return edge_has_vertex(e1, e2.first->vertex(e2.second))
&& edge_has_vertex(e1, e2.first->vertex(e2.third));
}
bool is_same_facet(const Facet& f1, const Facet& f2) const
{
return f1 == f2 || triangulation().mirror_facet(f1) == f2;
}
boost::optional<Vertex_handle> shared_vertex(const Edge& e1, const Edge& e2) const
{
Vertex_handle v1a = e1.first->vertex(e1.second);
Vertex_handle v1b = e1.first->vertex(e1.third);
@ -1065,22 +1013,22 @@ private:
return v1a;
else if (v1b == v2a || v1b == v2b)
return v1b;
std::cerr << "There is no vertex shared by e1 and e2" << std::endl;
CGAL_assertion(false);
return Vertex_handle();
else
return {};
}
Facet shared_facet(const Edge& e1, const Edge& e2) const
boost::optional<Facet> shared_facet(const Edge& e1, const Edge& e2) const
{
Vertex_handle v2a = e2.first->vertex(e2.second);
Vertex_handle v2b = e2.first->vertex(e2.third);
Vertex_handle sv = shared_vertex(e1, e2);
auto sv_opt = shared_vertex(e1, e2);
if(!sv_opt) return {};
Vertex_handle sv = *sv_opt;
Vertex_handle nsv2 = (sv == v2a) ? v2b : v2a;
typename Tr::Facet_circulator circ
= triangulation()->incident_facets(e1);
= triangulation().incident_facets(e1);
typename Tr::Facet_circulator end = circ;
do
{
@ -1092,33 +1040,30 @@ private:
}
} while (++circ != end);
std::cerr << "There is no facet shared by e1 and e2" << std::endl;
CGAL_assertion(false);
return Facet(Cell_handle(), 0);
return {};
}
Facet shared_facet(const Edge& e, const Vertex_handle v) const
{
typename Tr::Facet_circulator circ
= triangulation()->incident_facets(e);
= triangulation().incident_facets(e);
typename Tr::Facet_circulator end = circ;
do
{
Facet f = *circ;
int i;
if (triangulation()->has_vertex(f, v, i))
if (facet_has_vertex(f, v))
return f;
} while (++circ != end);
std::cerr << "There is no facet shared by e and v" << std::endl;
CGAL_assertion(false);
CGAL_unreachable();
return Facet(Cell_handle(), 0);
}
Cell_handle shared_cell(const Edge& e, const Vertex_handle v) const
{
typename Tr::Cell_circulator circ
= triangulation()->incident_cells(e);
= triangulation().incident_cells(e);
typename Tr::Cell_circulator end = circ;
do
{
@ -1128,7 +1073,7 @@ private:
} while (++circ != end);
std::cerr << "There is no cell shared by e and v" << std::endl;
CGAL_assertion(false);
CGAL_unreachable();
return Cell_handle();
}
@ -1145,6 +1090,11 @@ private:
}
}
Cell_handle shared_cell(const Edge e1, const Edge e2) const {
auto facet = shared_facet(e1, e2.first->vertex(e2.second));
return shared_cell(facet, e2.first->vertex(e2.third));
}
};//class Triangulation_segment_simplex_iterator_3
} // namespace CGAL

8
Triangulation_3/test/Triangulation_3/CMakeLists.txt
View File

@ -25,6 +25,14 @@ create_single_source_cgal_program("test_simplex_3.cpp")
create_single_source_cgal_program("test_simplex_iterator_3.cpp" )
create_single_source_cgal_program("test_segment_cell_traverser_3.cpp" )
create_single_source_cgal_program("test_segment_simplex_traverser_3.cpp" )
if(cxx_std_17 IN_LIST CMAKE_CXX_COMPILE_FEATURES)
target_compile_features(test_segment_simplex_traverser_3 PRIVATE cxx_std_17)
else()
message(
STATUS
"NOTICE: test_segment_simplex_traverser_3.cpp requires C++17 and will not be compiled."
)
endif()
create_single_source_cgal_program("test_static_filters.cpp")
create_single_source_cgal_program("test_triangulation_3.cpp")
create_single_source_cgal_program("test_io_triangulation_3.cpp")

452
Triangulation_3/test/Triangulation_3/test_segment_simplex_traverser_3.cpp
View File

@ -1,5 +1,8 @@
//#define CGAL_DEBUG_TRIANGULATION_SEGMENT_TRAVERSER_3 1
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Delaunay_triangulation_3.h>
#include <CGAL/Base_with_time_stamp.h>
#include <CGAL/IO/io.h>
#include <assert.h>
#include <iostream>
@ -15,18 +18,400 @@
// Define the kernel.
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point_3;
typedef Kernel::Segment_3 Segment_3;
typedef Kernel::Triangle_3 Triangle_3;
typedef Kernel::Tetrahedron_3 Tetrahedron_3;
// Define the structure.
typedef CGAL::Delaunay_triangulation_3<Kernel> DT;
typedef CGAL::Base_with_time_stamp<CGAL::Triangulation_vertex_base_3<Kernel>> Vb;
typedef CGAL::Delaunay_triangulation_cell_base_3<Kernel> Cb;
typedef CGAL::Triangulation_data_structure_3<Vb, Cb> Tds;
typedef CGAL::Delaunay_triangulation_3<Kernel, Tds> DT;
typedef DT::Cell_handle Cell_handle;
typedef DT::Edge Edge;
typedef DT::Facet Facet;
typedef DT::Vertex_handle Vertex_handle;
typedef DT::Simplex Simplex;
typedef DT::Segment_simplex_iterator Segment_simplex_iterator;
#include "test_triangulation_simplex_3_debug.h"
// a function to insert without spatial sorting
template <typename Point_it>
void insert(DT& dt, Point_it first, Point_it end) {
for(; first != end; ++first) {
dt.insert(*first);
}
}
static const std::vector<Point_3> bbox_points =
{
{ -10.1, -10, -10.08 },
{ -10.2, 10, -10.07 },
{ 10.3, 10, -10.06 },
{ 10.4, -10, -10.05 },
{ -10.5, -10, 10.04 },
{ -10.6, 10, 10.03 },
{ 10.7, 10, 10.02 },
{ 10.8, -10, 10.01 },
};
DT dt;
std::string result_string;
bool reverse_sort_vertex_handles(Vertex_handle v1, Vertex_handle v2) {
return v1.operator->() > v2.operator->();
};
auto vertices_of_simplex(Simplex simplex) -> std::array<Vertex_handle, 4> {
std::array<Vertex_handle, 4> vertices = { Vertex_handle{}, Vertex_handle{}, Vertex_handle{}, Vertex_handle{} };
switch(simplex.dimension()) {
case 0: {
vertices[0] = static_cast<Vertex_handle>(simplex);
break;
}
case 1: {
const auto [c, index1, index2] = static_cast<Edge>(simplex);
vertices[0] = c->vertex(index1);
vertices[1] = c->vertex(index2);
break;
}
case 2: {
const auto [c, index] = static_cast<Facet>(simplex);
vertices[0] = c->vertex(DT::vertex_triple_index(index, 0));
vertices[1] = c->vertex(DT::vertex_triple_index(index, 1));
vertices[2] = c->vertex(DT::vertex_triple_index(index, 2));
break;
}
case 3: {
const auto c = static_cast<Cell_handle>(simplex);
vertices[0] = c->vertex(0);
vertices[1] = c->vertex(1);
vertices[2] = c->vertex(2);
vertices[3] = c->vertex(3);
break;
}
default: CGAL_unreachable();
}
std::sort(vertices.begin(), vertices.end(), reverse_sort_vertex_handles);
for(int i = 0; i < 4; ++i) {
assert((i <= simplex.dimension()) == (vertices[i] != Vertex_handle{}));
}
return vertices;
}
std::variant<Point_3, Segment_3, Triangle_3, Tetrahedron_3> get_simplex_geometry(Simplex simplex) {
switch(simplex.dimension()) {
case 0: {
return static_cast<Vertex_handle>(simplex)->point();
}
case 1: {
const auto [c, index1, index2] = static_cast<Edge>(simplex);
return Segment_3(c->vertex(index1)->point(), c->vertex(index2)->point());
}
case 2: {
const auto [c, index] = static_cast<Facet>(simplex);
return Triangle_3(c->vertex(DT::vertex_triple_index(index, 0))->point(),
c->vertex(DT::vertex_triple_index(index, 1))->point(),
c->vertex(DT::vertex_triple_index(index, 2))->point());
}
case 3: {
const auto c = static_cast<Cell_handle>(simplex);
return Tetrahedron_3(c->vertex(0)->point(),
c->vertex(1)->point(),
c->vertex(2)->point(),
c->vertex(3)->point());
}
default: CGAL_unreachable();
}
}
void visit_simplex(Point_3 a, Point_3 b, Simplex s, std::optional<Simplex> previous_simplex_optional) {
auto d = s.dimension();
if(3 == d && dt.is_infinite(static_cast<Cell_handle>(s))) {
result_string += 'I';
} else {
result_string += std::to_string(d);
}
std::clog << debug_simplex(s) << '\n';
const bool does_intersect_ab = (3 == d && dt.is_infinite(static_cast<Cell_handle>(s))) ||
std::visit(
[&](auto geometry) { return CGAL::do_intersect(Segment_3(a, b), geometry); },
get_simplex_geometry(s));
if(!does_intersect_ab) {
CGAL_error_msg("the simplex does not intersect the query segment");
}
if(previous_simplex_optional) {
// this block checks that consecutive simplices are incident
using Set = std::array<Vertex_handle, 4>;
Set prev_vertices = vertices_of_simplex(*previous_simplex_optional);
Set s_vertices = vertices_of_simplex(s);
if(previous_simplex_optional->dimension() < s.dimension()) {
std::swap(prev_vertices, s_vertices);
std::swap(*previous_simplex_optional, s);
}
if(!std::includes(
prev_vertices.begin(), prev_vertices.begin() + 1 + previous_simplex_optional->dimension(),
s_vertices.begin(), s_vertices.begin() + 1 + s.dimension(),
reverse_sort_vertex_handles))
{
CGAL_error_msg("consecutive simplices are not incident");
}
}
};
bool test_vfefv(bool with_bbox = false)
{
std::cerr << "## test_vfefv(" << std::boolalpha << with_bbox << ")\n";
result_string.clear();
static const std::vector<Point_3> points =
{
{ -1, 0, 0 },
{ 0, 1, 0 },
{ 0, -1, 0 },
{ 5, 0, 0 },
{ 6, 2, 2 },
{ 6, -2, -2 },
};
dt.clear();
insert(dt, points.begin(), points.end());
if(with_bbox) insert(dt, bbox_points.begin(), bbox_points.end());
const auto v = dt.finite_vertex_handles().to<std::vector>();
Cell_handle c; int i, j, k;
assert(dt.is_facet(v[0], v[1], v[2], c, i, j, k));
assert(dt.is_facet(v[1], v[2], v[3], c, i, j, k));
assert(dt.is_cell (v[1], v[2], v[3], v[4], c));
assert(dt.is_cell (v[1], v[2], v[3], v[5], c));
std::optional<Simplex> previous{};
for(auto s: dt.segment_traverser_simplices(v[0], v[3])) {
visit_simplex(points[0], points[3], s, previous);
previous = s;
}
static const std::string expected_result_string = "02120";
bool ok = (result_string == expected_result_string);
if(!ok) {
std::cerr << "test_vfefv failed\n";
std::cerr << " result_string is " << result_string << " instead of "
<< expected_result_string << '\n';
}
return ok;
}
bool test_a_simple_tetrahedron(const std::vector<Point_3>& points) {
DT dt2;
for(const auto& p: points) dt2.insert(p);
bool ok = true;
auto test = [&](Point_3 a, Point_3 b, std::string expected_result) {
// This test function calls `do_test` with four configurations:
// - with [ab] and [ba],
// - and with or without a bbox around the central tetrahedron.
dt = dt2;
auto do_with_or_without_bbox = [&](Point_3 a, Point_3 b, bool with_bbox, std::string expected_result) {
auto do_it = [&](auto from, auto to) {
bool exception_thrown = false;
result_string.clear();
try {
std::optional<Simplex> previous_simplex;
for(auto s: dt.segment_traverser_simplices(from, to)) {
visit_simplex(a, b, s, previous_simplex);
previous_simplex = s;
}
} catch(const CGAL::Assertion_exception& e) {
CGAL::get_static_warning_handler()("Assertion", e.expression().c_str(),
e.filename().c_str(),
e.line_number(),
e.message().c_str());
exception_thrown = true;
}
if(result_string != expected_result || exception_thrown) {
std::cerr << "test_a_simple_tetrahedron failed on case " << expected_result
<< (with_bbox ? " with bbox\n" : "\n");
ok = false;
}
if(result_string != expected_result) {
std::cerr << " result_string is " << result_string << " instead of "
<< expected_result << '\n';
}
if(exception_thrown) {
std::cerr << " exception thrown\n";
}
}; // end do_it
std::clog << "### Case " << expected_result;
if(with_bbox) std::clog << " with bbox";
std::clog << '\n';
std::clog << "from (" << a << ") to (" << b << ")\n";
do_it(a, b);
// then re-test using vertex handles, if possible
Vertex_handle va{};
Vertex_handle vb{};
DT::Locate_type lt;
int i, j;
auto c = dt.locate(a, lt, i, j);
if(lt == DT::VERTEX) va = c->vertex(i);
c = dt.locate(b, lt, i, j);
if(lt == DT::VERTEX) vb = c->vertex(i);
if(va != Vertex_handle{} && vb != Vertex_handle{}) {
std::clog << "from vertex" << display_vert(va) << " to vertex" << display_vert(vb) << ")\n";
do_it(va, vb);
};
}; // end do_with_or_without_bbox
std::string expected_result_reversed{expected_result.rbegin(), expected_result.rend()};
do_with_or_without_bbox(a, b, false, expected_result);
do_with_or_without_bbox(b, a, false, expected_result_reversed);
std::replace(expected_result.begin(), expected_result.end(), 'I', '3');
std::replace(expected_result_reversed.begin(), expected_result_reversed.end(), 'I', '3');
insert(dt, bbox_points.begin(), bbox_points.end());
do_with_or_without_bbox(a, b, true, expected_result);
do_with_or_without_bbox(b, a, true, expected_result_reversed);
}; // end test() lambda
// [010] queries entering by a vertex and exiting by the other vertex of an incident edge,
// on the line (x,0,0)
test({ 0, 0, 0}, {.5, 0, 0}, "01");
test({ 0, 0, 0}, { 1, 0, 0}, "010");
test({ 0, 0, 0}, { 2, 0, 0}, "010I");
test({-1, 0, 0}, { 2, 0, 0}, "I010I");
test({-1, 0, 0}, { 1, 0, 0}, "I010");
test({-1, 0, 0}, {.5, 0, 0}, "I01");
// x [020] is not possible, because that would have passed through the edge -> see [010]
// [021] queries entering by a vertex and exiting by the opposite edge of a same face,
// on the line (x,x,0) (y==x)
test({ 0, 0, 0}, {.2, .2, 0}, "02");
test({ 0, 0, 0}, {.5, .5, 0}, "021");
test({ 0, 0, 0}, { 1, 1, 0}, "021I");
test({-1, -1, 0}, { 1, 1, 0}, "I021I");
test({-1, -1, 0}, {.5, .5, 0}, "I021");
test({-1, -1, 0}, {.2, .2, 0}, "I02");
// x [030] (entering by a vertex and exiting by a non-adjacent vertex) is not possible
// because that would have passed by the common edge -> see [010]
// x [031] (entering by a vertex and exiting by a non-incident edge), is not possible
// because that would have passed by the face -> see [021]
// [032] queries entering by a vertex and exiting by the opposite facet,
// on the line x==y==0.25-0.25z
test({ 0, 0, 1}, { .25, .25, .25 }, "03");
test({ 0, 0, 1}, { .25, .25, 0 }, "032");
test({ 0, 0, 1}, { .5, .5, -.1 }, "032I");
test({-.25,-.25, 2}, { .28125, .28125, -.125}, "I032I");
test({-.25,-.25, 2}, { .25, .25, 0 }, "I032");
test({-.25,-.25, 2}, { .125, .125, .5 }, "I03");
// [121] queries entering by an edge and exiting by an edge of the same face,
// on the line (x,.5,0)
test({0, .5, 0}, {.2, .5, 0}, "12");
test({0, .5, 0}, {.5 , .5, 0}, "121");
test({0, .5, 0}, {.6 , .5, 0}, "121I");
test({-.1, .5, 0}, {.6 , .5, 0}, "I121I");
test({-.1, .5, 0}, {.5 , .5, 0}, "I121");
test({-.1, .5, 0}, {.25, .5, 0}, "I12");
// x [130] (entering by an edge and exiting by a non-incident vertex) is not possible
// because that would have passed by the common face -> see [120] ([021] in reverse)
// [131] entering by an edge and exiting by the opposite edge in the cell,
// on the line x==y==0.5-z, also known as (.5-z, .5-z, z)
test({ 0, 0, .5 }, { .25, .25, .25 }, "13");
test({ 0, 0, .5 }, { .5, .5, 0 }, "131");
test({ 0, 0, .5 }, { .625, .625, -.125}, "131I");
test({ -.125, -.125, .625}, { .625, .625, -.125}, "I131I");
test({ -.125, -.125, .625}, { .5 , .5 , 0 }, "I131");
test({ -.125, -.125, .625}, { .25, .25, .25 }, "I13");
// [132] queries entering by an edge and exiting by a facet, on the line (x, .25-x, x)
test({ 0, .25, 0}, { .20, .05, .20}, "13");
test({ 0, .25, 0}, { .25, 0, .25}, "132");
test({ 0, .25, 0}, { .5 ,-.25, .5 }, "132I");
test({-.5, .75,-.5}, { .5 ,-.25, .5 }, "I132I");
test({-.5, .75,-.5}, { .25, 0, .25}, "I132");
test({-.5, .75,-.5}, { .20, .05, .20}, "I13");
// [232] queries entering by a facet and exiting by a facet, on the line (x,.5-x,.2)
test({ 0, .5, .2}, {.2, .3, .2}, "23");
test({ 0, .5, .2}, {.5, 0, .2}, "232");
test({ 0, .5, .2}, {.55, -.05, .2}, "232I");
test({-.05, .45, .2}, {.55, -.05, .2}, "I232I");
test({-.05, .45, .2}, {.5, 0, .2}, "I232");
test({-.05, .45, .2}, {.2, .3, .2}, "I23");
// special case: queries stay in a single simplex
test({ -.125, -.125, .625}, { -.125, -.125, .6251}, "I");
test({ .25, .25, .25}, { .20, .25, .25}, "3");
test({ 0, .5, .2}, { 0, .4, .2}, "2");
test({ 0, .5, .0}, { 0, .6, .0}, "1");
return ok;
}
bool test_a_simple_tetrahedron() {
bool ok = true;
std::cout << "## test_a_simple_tetrahedron()\n"
<< R"#(
This test uses with a trivial tetrahedron, and is launched with all the
24 permutations of the four vertices. There are 7 test cases, and for each of
them 6 different segments are tested:
- The segment starts at the incoming boundary of the tetrahedron and ends
inside.
- The segment starts at the incoming boundary of the tetrahedron and ends
at the outgoing boundary.
- The segment starts at the incoming boundary of the tetrahedron, goes out
and beyond.
- The segment starts before the tetrahedron, goes through it, and comes out.
- The segment starts before the tetrahedron and ends at the outgoing boundary.
- The segment starts before the tetrahedron and ends inside.
For each of those query segments, 4 tests are performed:
- with/without 8 extra vertices in the triangulation, forming a bounding
box,
- and in the direct and reverse direction.
In total, 4032 tests are performed...
)#";
std::cout.flush();
std::vector<Point_3> points = {
{0., 0., 0.},
{1., 0., 0.},
{0., 1., 0.},
{0., 0., 1.},
};
std::sort(points.begin(), points.end());
do {
std::cout << "### new permutation of the four points\n";
for(const auto& p: points) std::cout << " " << p << '\n';
std::cout << std::endl;
ok = test_a_simple_tetrahedron(points) && ok;
}
while (std::next_permutation(points.begin(), points.end()));
return ok;
}
bool test(const DT& dt,
const std::pair<Point_3, Point_3>& query,
const std::array<unsigned, 4>& expected_result);
int main(int, char* [])
{
std::cerr.precision(17);
std::clog.precision(17);
std::cout.precision(17);
bool ok = true;
ok = test_a_simple_tetrahedron() && ok;
const std::vector<Point_3> points = { { -2, 0, 0 },
{ 2, 0, 0 },
{ 0, 1, -1 },
@ -43,15 +428,13 @@ int main(int, char* [])
};
std::vector<DT::Vertex_handle> vertices;
vertices.reserve(points.size());
DT dt;
dt.clear();
for(auto p: points) vertices.push_back(dt.insert(p));
Cell_handle c;
assert(dt.is_valid());
assert(dt.is_cell(vertices[0], vertices[2], vertices[3], vertices[4], c));
assert(dt.is_cell(vertices[1], vertices[2], vertices[3], vertices[4], c));
std::cerr << dt.number_of_finite_cells() << '\n';
const std::vector < std::pair<Point_3, Point_3>> queries = {
{{-1, 0, 0}, { 1, 0, 0}}, // CFC
{{-1, 0, 0}, { 2, 0, 0}}, // CFCV
@ -74,12 +457,20 @@ int main(int, char* [])
{2, 1, 1, 0} // reverse case: FVEV
};
bool ok = true;
for(std::size_t i=0; i<queries.size(); ++i) {
if(!test(dt, queries[i], expected_results[i])) ok = false;
ok = test(dt, queries[i], expected_results[i]) && ok;
}
std::cout << "Done (" << queries.size() << " queries)\n";
return ok ? EXIT_SUCCESS : EXIT_FAILURE;
ok = test_vfefv() && ok;
ok = test_vfefv(true) && ok;
if(ok) {
std::cout << "All tests passed\n";
return EXIT_SUCCESS;
} else {
std::cout << "Some tests failed\n";
return EXIT_FAILURE;
}
}
bool test(const DT& dt,
@ -100,6 +491,7 @@ bool test(const DT& dt,
// Count the number of finite cells traversed.
unsigned int inf = 0, fin = 0;
std::optional<Simplex> previous;
for (; st != stend; ++st)
{
if (st->dimension() == 3
@ -108,46 +500,18 @@ bool test(const DT& dt,
else {
++fin;
visit_simplex(p1, p2, *st, previous);
previous = *st;
switch (st->dimension()) {
case 2: {
++nb_facets;
std::cout << "facet " << std::endl;
DT::Facet f = *st;
std::cout << " ( " << f.first->vertex((f.second + 1) & 3)->point()
<< " " << f.first->vertex((f.second + 2) & 3)->point()
<< " " << f.first->vertex((f.second + 3) & 3)->point()
<< " )\n";
break;
case 2: ++nb_facets; break;
case 1: ++nb_edges; break;
case 0: ++nb_vertex; break;
case 3: ++nb_cells; break;
default: CGAL_unreachable();
}
case 1: {
++nb_edges;
std::cout << "edge " << std::endl;
DT::Edge e = *st;
std::cout << " ( " << e.first->vertex(e.second)->point() << " "
<< e.first->vertex(e.third)->point() << " )\n";
break;
}
case 0: {
++nb_vertex;
std::cout << "vertex " << std::endl;
DT::Vertex_handle v = *st;
std::cout << " ( " << v->point() << " )\n";
break;
}
case 3: {
++nb_cells;
std::cout << "cell \n ( ";
DT::Cell_handle ch = *st;
for (int i = 0; i < 4; ++i)
std::cout << ch->vertex(i)->point() << " ";
std::cout << " )\n";
break;
}
default:
CGAL_unreachable();
} // end switch
}
}
}
#ifdef CGAL_TRIANGULATION_3_VERBOSE_TRAVERSER_EXAMPLE
std::cout << "While traversing from " << st.source()

14
Triangulation_3/test/Triangulation_3/test_simplex_iterator_3.cpp
View File

@ -1,5 +1,7 @@
// #define CGAL_DEBUG_TRIANGULATION_SEGMENT_TRAVERSER_3 1
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Delaunay_triangulation_3.h>
#include <CGAL/Base_with_time_stamp.h>
#include <assert.h>
#include <iostream>
@ -17,14 +19,20 @@ typedef Kernel::Vector_3 Vector_3;
typedef Kernel::Segment_3 Segment_3;
// Define the structure.
typedef CGAL::Delaunay_triangulation_3< Kernel > DT;
typedef CGAL::Base_with_time_stamp<CGAL::Triangulation_vertex_base_3<Kernel>> Vb;
typedef CGAL::Delaunay_triangulation_cell_base_3<Kernel> Cb;
typedef CGAL::Triangulation_data_structure_3<Vb, Cb> Tds;
typedef CGAL::Delaunay_triangulation_3< Kernel, Tds > DT;
typedef DT::Vertex_handle Vertex_handle;
typedef DT::Cell_handle Cell_handle;
typedef DT::Edge Edge;
typedef DT::Facet Facet;
typedef DT::Simplex Simplex;
typedef DT::Segment_simplex_iterator Segment_simplex_iterator;
#include "test_triangulation_simplex_3_debug.h"
void test_vertex_edge_vertex(const DT& dt, const std::size_t& nb_tests)
{
std::cout << "* test_vertex_edge_vertex *" << std::endl;
@ -299,8 +307,8 @@ void test_triangulation_on_a_grid()
unsigned int nb_facets = 0, nb_edges = 0, nb_vertex = 0;
for (; st != st.end(); ++st)
{
std::cout << st->dimension() << " ";
std::cout.flush();
std::cerr << st->dimension() << " ";
std::cerr << debug_simplex(*st) <<'\n';
if (st->dimension() == 3)
{
if (dt.is_infinite(Cell_handle(*st))) ++inf;

68
Triangulation_3/test/Triangulation_3/test_triangulation_simplex_3_debug.h Normal file
View File

@ -0,0 +1,68 @@
template <typename Vertex_handle>
auto display_vert(Vertex_handle v) {
std::stringstream os;
os.precision(17);
if(v->time_stamp() == 0) {
os << "inf";
} else {
os << '#' << v->time_stamp() << "=(" << v->point() << ")";
}
return os.str();
};
template <typename DT>
struct Debug_simplex {
using Cell_handle = typename DT::Cell_handle;
using Edge = typename DT::Edge;
using Facet = typename DT::Facet;
using Vertex_handle = typename DT::Vertex_handle;
using Simplex = typename DT::Simplex;
Simplex simplex;
template<typename Dt, typename CharT, typename Traits>
friend
std::basic_ostream<CharT, Traits>&
operator<<(std::basic_ostream<CharT, Traits>& os, const Debug_simplex<Dt>& d) {
auto&& simplex = d.simplex;
switch(simplex.dimension()) {
case 0: {
os << "- vertex " << display_vert(static_cast<Vertex_handle>(simplex));
break;
}
case 1: {
const auto [c, index1, index2] = static_cast<Edge>(simplex);
os << "- edge "
<< display_vert(c->vertex(index1)) << " - "
<< display_vert(c->vertex(index2));
break;
}
case 2: {
const auto [c, index] = static_cast<Facet>(simplex);
os << "- facet "
<< display_vert(c->vertex(DT::vertex_triple_index(index, 0))) << " - "
<< display_vert(c->vertex(DT::vertex_triple_index(index, 1))) << " - "
<< display_vert(c->vertex(DT::vertex_triple_index(index, 2)));
break;
}
case 3: {
const auto c = static_cast<Cell_handle>(simplex);
os << "- cell "
<< display_vert(c->vertex(0)) << " - "
<< display_vert(c->vertex(1)) << " - "
<< display_vert(c->vertex(2)) << " - "
<< display_vert(c->vertex(3));
break;
}
default: CGAL_assume(false);
}
return os;
};
};
#include <CGAL/Triangulation_simplex_3.h>
template <typename Triangulation>
auto debug_simplex(CGAL::Triangulation_simplex_3<Triangulation> simplex) {
return Debug_simplex<Triangulation>{simplex};
}