diff --git a/STL_Extension/include/CGAL/Iterator_range.h b/STL_Extension/include/CGAL/Iterator_range.h index e1d41818d5a7..e42441e98588 100644 --- a/STL_Extension/include/CGAL/Iterator_range.h +++ b/STL_Extension/include/CGAL/Iterator_range.h @@ -78,6 +78,13 @@ namespace CGAL { { return std::tuple{this->first, this->second}; } + + template class Container> + auto to() const + { + using V = std::remove_cv_t>; + return Container(begin(), end()); + } }; template diff --git a/TDS_3/include/CGAL/Triangulation_simplex_3.h b/TDS_3/include/CGAL/Triangulation_simplex_3.h index aa2109b89ced..0754fb021c59 100644 --- a/TDS_3/include/CGAL/Triangulation_simplex_3.h +++ b/TDS_3/include/CGAL/Triangulation_simplex_3.h @@ -96,8 +96,10 @@ class Triangulation_simplex_3 { // 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 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 s0, } return false; case (3): - return (&(*s0.ch) == &(*s1.ch)); + return s0.ch.operator->() == s1.ch.operator->(); } CGAL_error(); return false; diff --git a/Triangulation_3/include/CGAL/Triangulation_3/internal/Triangulation_segment_traverser_3_impl.h b/Triangulation_3/include/CGAL/Triangulation_3/internal/Triangulation_segment_traverser_3_impl.h index b93c9a3f8745..8040af90f2cd 100644 --- a/Triangulation_3/include/CGAL/Triangulation_3/internal/Triangulation_segment_traverser_3_impl.h +++ b/Triangulation_3/include/CGAL/Triangulation_3/internal/Triangulation_segment_traverser_3_impl.h @@ -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 Triangulation_segment_cell_iterator_3::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:: 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::Simplex, Triangulation_segment_cell_iterator_3::walk_to_next_3(const Simplex& prev, const Simplex& cur) const { - std::array 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()) }; - - int inside=0,outside=0,regular_case=0,degenerate=0; - Cell_handle nnext; - - 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; + const auto cur_cell = cur.cell; + std::array 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; + + 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; + Simplex cur_after_walk; + + 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 } - else { - if (_tr->orientation(*vert[i], *vert[j1], *vert[j2], _target) == POSITIVE) { - nnext = std::get<0>(cur)->neighbor(j0); - outside = j0; - } - else - inside = 2; + } // end o1 <= 0 + else + { // 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 } } - 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; + } // 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 { - if (_tr->orientation(*vert[i], *vert[j1], *vert[j2], _target) == POSITIVE) { - nnext = std::get<0>(cur)->neighbor(j0); - outside = j0; - } - else - inside = 4; + 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 } } - 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 + } // 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 - inside = 5; } - else { - if (_tr->orientation(*vert[i], *vert[j1], *vert[j2], _target) == POSITIVE) { - nnext = std::get<0>(cur)->neighbor(j0); - outside = j0; + } 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; } - 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; + if(!degenerate) { + return {prev_after_walk, cur_after_walk}; } + } - // 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; + // 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 o; + std::array 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; } - 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; - } + } + } 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; + } - 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]; - } + // For the remembering stochastic walk, we start trying with a random facet. + CGAL_assertion_code(bool incell = true;) - if( Or == 0 ) { - // Either the target is not inside the pyramid, - // or the pyramid is degenerate. - vert[li] = backup; + 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]; + } - // 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); + if(Or == 0) { + // Either the target is not inside the pyramid, + // or the pyramid is degenerate. + vert[li] = backup_vert_li; + continue; + } - 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 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); } - - // The target lies inside this cell. - 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; - - 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 ); - 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; + 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: - prev_after_walk = Simplex( std::get<0>(cur), Tr::OUTSIDE_AFFINE_HULL, -1, -1 ); - CGAL_assertion( false ); + CGAL_unreachable(); + return std::make_pair(prev, cur); } + CGAL_unreachable(); + } - Simplex cur_after_walk(Cell_handle(), Tr::OUTSIDE_AFFINE_HULL, -1, -1); - return std::make_pair(prev_after_walk, cur_after_walk); + // 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::opposite_edge( case 5: return Edge(c, 2, 3); } - CGAL_assertion(false); + CGAL_unreachable(); return Edge(); } diff --git a/Triangulation_3/include/CGAL/Triangulation_segment_traverser_3.h b/Triangulation_3/include/CGAL/Triangulation_segment_traverser_3.h index 904421786ba2..b7ccdeb6a0e0 100644 --- a/Triangulation_3/include/CGAL/Triangulation_segment_traverser_3.h +++ b/Triangulation_3/include/CGAL/Triangulation_segment_traverser_3.h @@ -17,7 +17,6 @@ #include #include -#include #include #include @@ -27,6 +26,7 @@ #include #include +#include // If defined, type casting is done statically, // reducing type-safety overhead. @@ -113,7 +113,13 @@ class Triangulation_segment_cell_iterator_3 typedef typename Tr::Locate_type Locate_type; //< defines the simplex type returned from location. - typedef std::tuple 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 @@ class Triangulation_segment_cell_iterator_3 template < class Tr2, class Inc2 > struct Rebind { typedef Triangulation_segment_cell_iterator_3 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 @@ class Triangulation_segment_cell_iterator_3 */ 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 @@ class Triangulation_segment_cell_iterator_3 */ Cell* operator->() { - return &*std::get<0>(_cur); + return &*(_cur.cell); } // provides an indirection operator. @@ -265,7 +322,7 @@ class Triangulation_segment_cell_iterator_3 */ Cell& operator*() { - return *std::get<0>(_cur); + return *(_cur.cell); } // provides a conversion operator. @@ -273,7 +330,7 @@ class Triangulation_segment_cell_iterator_3 */ operator Cell_handle() const { - return std::get<0>(_cur); + return _cur.cell; } // provides a conversion operator. @@ -299,6 +356,10 @@ class Triangulation_segment_cell_iterator_3 { lt = this->lt(); li = this->li(); lj = this->lj(); } + std::tuple 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 @@ class Triangulation_segment_cell_iterator_3 { lt = prev_lt(); li = prev_li(); lj = prev_lj(); } + std::tuple 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 @@ class Triangulation_segment_cell_iterator_3 */ 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 @@ class Triangulation_segment_cell_iterator_3 */ bool operator!=( const Cell_handle& ch ) const { - return ch != std::get<0>(_cur); + return ch != _cur.cell; } // \} @@ -439,32 +504,33 @@ class Triangulation_segment_cell_iterator_3 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 @@ class Triangulation_segment_simplex_iterator_3 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 @@ class Triangulation_segment_simplex_iterator_3 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 @@ class Triangulation_segment_simplex_iterator_3 _curr_simplex = cell; break; default: - CGAL_assertion(false); + CGAL_unreachable(); }; } @@ -637,282 +706,145 @@ class Triangulation_segment_simplex_iterator_3 // 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 = 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(); - - 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(); - else - _curr_simplex = ch; - break; - } - - case Locate_type::EDGE: - if (facet_has_edge(get_facet(), Edge(chnext, linext, ljnext))) - set_curr_simplex_to_entry(); - else - _curr_simplex = ch; - break; - case Locate_type::FACET: - _curr_simplex = ch; - break; + 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); - case Locate_type::OUTSIDE_AFFINE_HULL: - _curr_simplex = Simplex_3(); - break; - - default: - CGAL_assertion(false); - }; - break; + if(_curr_simplex.dimension() == 3) { + set_curr_simplex_to_entry(); + return *this; } - case 1:/*Edge*/ + 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 (ch == _cell_iterator.previous()) - { - _curr_simplex = Simplex_3(); - break; - } - 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; - } - } + 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 (edge_has_vertex(get_edge(), chnext->vertex(linext))) - _curr_simplex = chnext->vertex(linext); - else - _curr_simplex = shared_facet(get_edge(), chnext->vertex(linext)); - break; + switch(_curr_simplex.dimension()) { + case 2: { /*Facet*/ + CGAL_assertion((ch_next == Cell_handle()) == (_cell_iterator == _cell_iterator.end())); - 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; + 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 = 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); + _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 = Simplex_3(); //edge-outside - 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);//should not happen - }; - 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 + CGAL_unreachable(); } - 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; + } 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 - { - 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 { + auto facet_opt = shared_facet(get_edge(), e_prev); + if(static_cast(facet_opt)) { + _curr_simplex = *facet_opt; } - } - 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); + else { + _curr_simplex = shared_cell(get_edge(), e_prev); } } - 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; + } 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 = ch_prev; //query goes through the cell + } break; + default: + CGAL_unreachable(); } - - 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 + } break; + case 0 :/*Vertex_handle*/ + { + 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 @@ class Triangulation_segment_simplex_iterator_3 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 @@ class Triangulation_segment_simplex_iterator_3 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 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 @@ class Triangulation_segment_simplex_iterator_3 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 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 @@ class Triangulation_segment_simplex_iterator_3 } } 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 @@ class Triangulation_segment_simplex_iterator_3 } 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 @@ class Triangulation_segment_simplex_iterator_3 } } + 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 diff --git a/Triangulation_3/test/Triangulation_3/CMakeLists.txt b/Triangulation_3/test/Triangulation_3/CMakeLists.txt index ffd521512de0..31ecea0cef42 100644 --- a/Triangulation_3/test/Triangulation_3/CMakeLists.txt +++ b/Triangulation_3/test/Triangulation_3/CMakeLists.txt @@ -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") diff --git a/Triangulation_3/test/Triangulation_3/test_segment_simplex_traverser_3.cpp b/Triangulation_3/test/Triangulation_3/test_segment_simplex_traverser_3.cpp index eb1fa780be49..45e4f46c318d 100644 --- a/Triangulation_3/test/Triangulation_3/test_segment_simplex_traverser_3.cpp +++ b/Triangulation_3/test/Triangulation_3/test_segment_simplex_traverser_3.cpp @@ -1,5 +1,8 @@ +//#define CGAL_DEBUG_TRIANGULATION_SEGMENT_TRAVERSER_3 1 #include #include +#include +#include #include #include @@ -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 DT; +typedef CGAL::Base_with_time_stamp> Vb; +typedef CGAL::Delaunay_triangulation_cell_base_3 Cb; +typedef CGAL::Triangulation_data_structure_3 Tds; +typedef CGAL::Delaunay_triangulation_3 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 +void insert(DT& dt, Point_it first, Point_it end) { + for(; first != end; ++first) { + dt.insert(*first); + } +} + +static const std::vector 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 { + std::array vertices = { Vertex_handle{}, Vertex_handle{}, Vertex_handle{}, Vertex_handle{} }; + switch(simplex.dimension()) { + case 0: { + vertices[0] = static_cast(simplex); + break; + } + case 1: { + const auto [c, index1, index2] = static_cast(simplex); + vertices[0] = c->vertex(index1); + vertices[1] = c->vertex(index2); + break; + } + case 2: { + const auto [c, index] = static_cast(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(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 get_simplex_geometry(Simplex simplex) { + switch(simplex.dimension()) { + case 0: { + return static_cast(simplex)->point(); + } + case 1: { + const auto [c, index1, index2] = static_cast(simplex); + return Segment_3(c->vertex(index1)->point(), c->vertex(index2)->point()); + } + case 2: { + const auto [c, index] = static_cast(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(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 previous_simplex_optional) { + auto d = s.dimension(); + if(3 == d && dt.is_infinite(static_cast(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(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; + 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 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(); + + 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 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& 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 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 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& query, const std::array& 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 points = { { -2, 0, 0 }, { 2, 0, 0 }, { 0, 1, -1 }, @@ -43,15 +428,13 @@ int main(int, char* []) }; std::vector 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> 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 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() diff --git a/Triangulation_3/test/Triangulation_3/test_simplex_iterator_3.cpp b/Triangulation_3/test/Triangulation_3/test_simplex_iterator_3.cpp index 048581c0d906..506fd81022cf 100644 --- a/Triangulation_3/test/Triangulation_3/test_simplex_iterator_3.cpp +++ b/Triangulation_3/test/Triangulation_3/test_simplex_iterator_3.cpp @@ -1,5 +1,7 @@ +// #define CGAL_DEBUG_TRIANGULATION_SEGMENT_TRAVERSER_3 1 #include #include +#include #include #include @@ -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> Vb; +typedef CGAL::Delaunay_triangulation_cell_base_3 Cb; +typedef CGAL::Triangulation_data_structure_3 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; diff --git a/Triangulation_3/test/Triangulation_3/test_triangulation_simplex_3_debug.h b/Triangulation_3/test/Triangulation_3/test_triangulation_simplex_3_debug.h new file mode 100644 index 000000000000..5ea79279b1cb --- /dev/null +++ b/Triangulation_3/test/Triangulation_3/test_triangulation_simplex_3_debug.h @@ -0,0 +1,68 @@ +template +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 +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 + friend + std::basic_ostream& + operator<<(std::basic_ostream& os, const Debug_simplex
& d) { + auto&& simplex = d.simplex; + switch(simplex.dimension()) { + case 0: { + os << "- vertex " << display_vert(static_cast(simplex)); + break; + } + case 1: { + const auto [c, index1, index2] = static_cast(simplex); + os << "- edge " + << display_vert(c->vertex(index1)) << " - " + << display_vert(c->vertex(index2)); + break; + } + case 2: { + const auto [c, index] = static_cast(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(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 + +template +auto debug_simplex(CGAL::Triangulation_simplex_3 simplex) { + return Debug_simplex{simplex}; +}