DGtal 1.3.0
Loading...
Searching...
No Matches
HalfEdgeDataStructure.h
1
17#pragma once
18
31#if defined(HalfEdgeDataStructure_RECURSES)
32#error Recursive header files inclusion detected in HalfEdgeDataStructure.h
33#else // defined(HalfEdgeDataStructure_RECURSES)
35#define HalfEdgeDataStructure_RECURSES
36
37#if !defined HalfEdgeDataStructure_h
39#define HalfEdgeDataStructure_h
40
42// Inclusions
43#include <iostream>
44#include <array>
45#include "DGtal/base/Common.h"
47
48namespace DGtal
49{
50
57 static std::size_t const HALF_EDGE_INVALID_INDEX = boost::integer_traits<std::size_t>::const_max;
58
60 // class HalfEdgeDataStructure
83 {
84 public:
85
87 typedef std::size_t Size;
89 typedef std::size_t Index;
98
99 // Defines the invalid index.
100 // JOL: works with Apple LLVM version 8.1.0 (clang-802.0.42), but does not work with gcc 4.8.4
101 //
102 // static constexpr Index const& INVALID_INDEX = boost::integer_traits<Index>::const_max;
103 //
104 // JOL: does NOT work with Apple LLVM version 8.1.0 (clang-802.0.42)
105 // (undefined reference at link time).
106 // static constexpr Index const INVALID_INDEX = boost::integer_traits<Index>::const_max;
107 // static Index const INVALID_INDEX = boost::integer_traits<Index>::const_max;
108 // BOOST_STATIC_CONSTANT( Index, INVALID_INDEX = boost::integer_traits<Index>::const_max );
109
110 typedef std::vector<HalfEdgeIndex> HalfEdgeIndexRange;
111 typedef std::vector<VertexIndex> VertexIndexRange;
112 typedef std::vector<EdgeIndex> EdgeIndexRange;
113 typedef std::vector<FaceIndex> FaceIndexRange;
114
116 typedef std::pair<VertexIndex, VertexIndex> Arc;
117 // A map from an arc (a std::pair of VertexIndex's) to its
118 // half edge index (i.e. and offset into the 'halfedge' sequence).
119 typedef std::map< Arc, Index > Arc2Index;
120 // A map from an arc (a std::pair of VertexIndex's) to its face
121 // index.
122 typedef std::map< Arc, FaceIndex > Arc2FaceIndex;
123
126 struct Edge
127 {
130
131 VertexIndex& start() { return v[0]; }
132 const VertexIndex& start() const { return v[0]; }
133
134 VertexIndex& end() { return v[1]; }
135 const VertexIndex& end() const { return v[1]; }
136
138 {
139 v[0] = v[1] = -1;
140 }
142 {
143 if ( vi <= vj ) { v[0] = vi; v[1] = vj; }
144 else { v[0] = vj; v[1] = vi; }
145 }
146 bool operator<( const Edge& other ) const
147 {
148 return ( start() < other.start() )
149 || ( ( start() == other.start() ) && ( end() < other.end() ) );
150 }
151 };
152
154 struct Triangle
155 {
157 std::array<VertexIndex,3> v;
158
159 VertexIndex& i() { return v[0]; }
160 const VertexIndex& i() const { return v[0]; }
161
162 VertexIndex& j() { return v[1]; }
163 const VertexIndex& j() const { return v[1]; }
164
165 VertexIndex& k() { return v[2]; }
166 const VertexIndex& k() const { return v[2]; }
167
169 {
170 v[0] = v[1] = v[2] = -1;
171 }
173 {
174 v[0] = v0;
175 v[1] = v1;
176 v[2] = v2;
177 }
178 };
179
182 typedef std::vector<VertexIndex> PolygonalFace;
183
188 struct HalfEdge
189 {
200
208 {}
209 };
210
211 public:
214
231 ( const std::vector<Triangle>& triangles, std::vector< Edge >& edges_out )
232 {
233 typedef std::set< Edge > EdgeSet;
234 typedef std::set< VertexIndex > VertexIndexSet;
235 VertexIndexSet vertexSet;
236 EdgeSet edgeSet;
237 for( const Triangle& T : triangles )
238 {
239 edgeSet.insert( Edge( T.i(), T.j() ) );
240 edgeSet.insert( Edge( T.j(), T.k() ) );
241 edgeSet.insert( Edge( T.k(), T.i() ) );
242 vertexSet.insert( T.i() );
243 vertexSet.insert( T.j() );
244 vertexSet.insert( T.k() );
245 }
246 edges_out.resize( edgeSet.size() );
247 Size e = 0;
248 for ( const Edge& edge : edgeSet )
249 {
250 edges_out.at(e) = edge;
251 ++e;
252 }
253 return vertexSet.size();
254 }
255
272 ( const std::vector<PolygonalFace>& polygonal_faces, std::vector< Edge >& edges_out );
273
296 bool build( const Size num_vertices,
297 const std::vector<Triangle>& triangles,
298 const std::vector<Edge>& edges );
299
322 bool build( const Size num_vertices,
323 const std::vector<PolygonalFace>& polygonal_faces,
324 const std::vector<Edge>& edges );
325
334 bool build( const std::vector<Triangle>& triangles )
335 {
336 std::vector<Edge> edges;
337 const Size nbVtx = getUnorderedEdgesFromTriangles( triangles, edges );
338 return build( nbVtx, triangles, edges );
339 }
340
349 bool build( const std::vector<PolygonalFace>& polygonal_faces )
350 {
351 std::vector<Edge> edges;
352 const Size nbVtx = getUnorderedEdgesFromPolygonalFaces( polygonal_faces, edges );
353 return build( nbVtx, polygonal_faces, edges );
354 }
355
357 void clear()
358 {
359 myHalfEdges.clear();
360 myVertexHalfEdges.clear();
361 myFaceHalfEdges.clear();
362 myEdgeHalfEdges.clear();
363 myArc2Index.clear();
364 }
365
367 Size nbHalfEdges() const { return myHalfEdges.size(); }
368
370 Size nbVertices() const { return myVertexHalfEdges.size(); }
371
373 Size nbEdges() const { return myEdgeHalfEdges.size(); }
374
376 Size nbFaces() const { return myFaceHalfEdges.size(); }
377
379 long Euler() const
380 { return (long) nbVertices() - (long) nbEdges() + (long) nbFaces(); }
381
384 const HalfEdge& halfEdge( const Index i ) const { return myHalfEdges.at( i ); }
385
389 {
390 const HalfEdge& he = myHalfEdges[ i ];
391 return std::make_pair( myHalfEdges[ he.opposite ].toVertex, he.toVertex );
392 }
393
398 { return halfEdgeIndexFromArc( std::make_pair( i, j ) ); }
399
402 Index halfEdgeIndexFromArc( const Arc& arc ) const
403 {
404 auto result = myArc2Index.find( arc );
405 return ( result == myArc2Index.end() ) ? HALF_EDGE_INVALID_INDEX : result->second;
406 }
407
419 { return findHalfEdgeIndexFromArc( std::make_pair( i, j ) ); }
420
431 {
432 const Index s = halfEdgeIndexFromVertexIndex( arc.first );
433 Index i = s;
434 do {
435 const HalfEdge& he = myHalfEdges[ i ];
436 if ( he.toVertex == arc.second ) return i;
437 i = halfEdge( he.opposite ).next;
438 } while ( i != s );
440 }
441
445 { return myVertexHalfEdges[ vi ]; }
446
450 { return myFaceHalfEdges[ fi ]; }
451
455 { return myEdgeHalfEdges[ ei ]; }
456
460 {
461 result.clear();
462 const Index start_hei = halfEdgeIndexFromVertexIndex( vi );
463 Index hei = start_hei;
464 do
465 {
466 const HalfEdge& he = halfEdge( hei );
467 result.push_back( he.toVertex );
468 hei = halfEdge( he.opposite ).next;
469 }
470 while ( hei != start_hei );
471 }
472
476 {
477 VertexIndexRange result;
478 getNeighboringVertices( vi, result );
479 return result;
480 }
481
485 {
486 Size nb = 0;
487 const Index start_hei = halfEdgeIndexFromVertexIndex( vi );
488 Index hei = start_hei;
489 do
490 {
491 const HalfEdge& he = halfEdge( hei );
492 nb++;
493 hei = halfEdge( he.opposite ).next;
494 }
495 while ( hei != start_hei );
496 return nb;
497 }
498
501 void getNeighboringFaces( const VertexIndex vi, FaceIndexRange& result ) const
502 {
503 result.clear();
504 const Index start_hei = halfEdgeIndexFromVertexIndex( vi );
505 Index hei = start_hei;
506 do
507 {
508 const HalfEdge& he = halfEdge( hei );
509 if( HALF_EDGE_INVALID_INDEX != he.face ) result.push_back( he.face );
510 hei = halfEdge( he.opposite ).next;
511 }
512 while ( hei != start_hei );
513 }
514
518 {
519 FaceIndexRange result;
520 getNeighboringFaces( vi, result );
521 return result;
522 }
523
526 bool isVertexBoundary( const VertexIndex vi ) const
527 {
529 }
530
537 {
538 VertexIndexRange result;
539 // std::set< VertexIndex > result;
540 for( Index hei = 0; hei < myHalfEdges.size(); ++hei )
541 {
542 const HalfEdge& he = halfEdge( hei );
543 if( HALF_EDGE_INVALID_INDEX == he.face )
544 result.push_back( he.toVertex );
545 }
546 return result;
547 }
548
552 std::vector< Index > boundaryHalfEdgeIndices() const
553 {
554 std::vector< Index > result;
555 for( Index hei = 0; hei < myHalfEdges.size(); ++hei )
556 {
557 const HalfEdge& he = halfEdge( hei );
558 if( HALF_EDGE_INVALID_INDEX == he.face )
559 result.push_back( hei );
560 }
561 return result;
562 }
566 std::vector< Arc > boundaryArcs() const
567 {
568 std::vector< Arc > result;
569 for( Index hei = 0; hei < myHalfEdges.size(); ++hei )
570 {
571 const HalfEdge& he = halfEdge( hei );
572 if( HALF_EDGE_INVALID_INDEX == he.face )
573 result.push_back( arcFromHalfEdgeIndex( hei ) );
574 }
575 return result;
576 }
577
580 Size nbSides( Index hei ) const
581 {
582 ASSERT( hei != HALF_EDGE_INVALID_INDEX );
583 Size nb = 0;
584 const Index start = hei;
585 do {
586 hei = halfEdge( hei ).next;
587 nb++;
588 }
589 while ( hei != start );
590 return nb;
591 }
592
595 Size nbSidesOfFace( const FaceIndex f ) const
596 {
597 return nbSides( halfEdgeIndexFromFaceIndex( f ) );
598 }
599
603 {
604 VertexIndexRange result;
605 result.reserve( 3 );
607 const Index start = hei;
608 do {
609 const HalfEdge& he = halfEdge( hei );
610 result.push_back( he.toVertex );
611 hei = he.next;
612 }
613 while ( hei != start );
614 return result;
615 }
616
617 // -------------------- Modification services -------------------------
618 public:
619
628 bool isFlippable( const Index hei ) const
629 {
630 ASSERT( hei != HALF_EDGE_INVALID_INDEX );
631 const HalfEdge& he = halfEdge( hei );
632 const Index hei2 = he.opposite;
633 const HalfEdge& he2 = halfEdge( hei2 );
634 // check if hei borders an infinite face.
635 if ( he.face == HALF_EDGE_INVALID_INDEX ) return false;
636 // check if hei2 borders an infinite face.
637 if ( he2.face == HALF_EDGE_INVALID_INDEX ) return false;
638 // Checks that he1 and he2 border a triangle.
639 if ( ( nbSides( hei ) != 3 ) || ( nbSides( hei2 ) != 3 ) ) return false;
640 // Checks that the two other vertices of the two surrounding
641 // triangles are not already neighbors
642 const VertexIndex v1 = halfEdge( he.next ).toVertex;
643 const VertexIndex v2 = halfEdge( he2.next ).toVertex;
644 const auto neighb_v1 = neighboringVertices( v1 );
645 const auto it_v2 = std::find( neighb_v1.cbegin(), neighb_v1.cend(), v2 );
646 return it_v2 == neighb_v1.cend();
647 }
648
664 void flip( const Index hei, bool update_arc2index = true )
665 {
666 const Index i1 = hei;
667 HalfEdge& he1 = myHalfEdges[ i1 ];
668 const Index i2 = he1.opposite;
669 HalfEdge& he2 = myHalfEdges[ i2 ];
670 const VertexIndex v2 = he1.toVertex;
671 const VertexIndex v1 = he2.toVertex;
672 const Index i1_next = he1.next;
673 const Index i2_next = he2.next;
674 HalfEdge& he1_next = myHalfEdges[ i1_next ];
675 HalfEdge& he2_next = myHalfEdges[ i2_next ];
676 const Index i1_next2 = he1_next.next;
677 const Index i2_next2 = he2_next.next;
678 myHalfEdges[ i1_next2 ].next = i2_next;
679 myHalfEdges[ i2_next2 ].next = i1_next;
680 he2_next.next = i1;
681 he1_next.next = i2;
682 he2_next.face = he1.face;
683 he1_next.face = he2.face;
684 he1.next = i1_next2;
685 he2.next = i2_next2;
686 he1.toVertex = he1_next.toVertex;
687 he2.toVertex = he2_next.toVertex;
688 // Reassign the mapping vertex v -> index of half edge starting from v
689 // (JOL): must check before reassign for boundary vertices.
690 if ( myVertexHalfEdges[ v1 ] == i1 ) myVertexHalfEdges[ v1 ] = i2_next;
691 if ( myVertexHalfEdges[ v2 ] == i2 ) myVertexHalfEdges[ v2 ] = i1_next;
692 // Reassign the mapping face f -> index of half edge contained in f
693 myFaceHalfEdges[ he1.face ] = i1;
694 myFaceHalfEdges[ he2.face ] = i2;
695 // No need to reassign edge... it has just changed of vertices
696 // but is still based on half-edges i1 and i2
697 // Now taking care of mapping (vertex,vertex) -> half-edge.
698 if ( update_arc2index )
699 {
700 myArc2Index.erase( Arc( v1, v2 ) );
701 myArc2Index.erase( Arc( v2, v1 ) );
702 const VertexIndex v2p = he1.toVertex;
703 const VertexIndex v1p = he2.toVertex;
704 myArc2Index[ Arc( v1p, v2p ) ] = i1;
705 myArc2Index[ Arc( v2p, v1p ) ] = i2;
706 }
707 }
708
727 VertexIndex split( const Index i, bool update_arc2index = true )
728 {
729 Index new_hei = myHalfEdges.size();
730 VertexIndex new_vtxi = myVertexHalfEdges.size();
731 EdgeIndex new_edgei = myEdgeHalfEdges.size();
732 FaceIndex new_facei = myFaceHalfEdges.size();
733 myHalfEdges.resize( new_hei + 6 );
734 myVertexHalfEdges.push_back( new_hei );
735 HalfEdge& hei = myHalfEdges[ i ];
736 HalfEdge& hei_next = myHalfEdges[ hei.next ];
737 Index j = hei.opposite;
738 HalfEdge& hej = myHalfEdges[ j ];
739 HalfEdge& hej_next = myHalfEdges[ hej.next ];
740 HalfEdge& he0 = myHalfEdges[ new_hei ];
741 HalfEdge& he1 = myHalfEdges[ new_hei + 1 ];
742 HalfEdge& he2 = myHalfEdges[ new_hei + 2 ];
743 HalfEdge& he3 = myHalfEdges[ new_hei + 3 ];
744 HalfEdge& he4 = myHalfEdges[ new_hei + 4 ];
745 HalfEdge& he5 = myHalfEdges[ new_hei + 5 ];
746 // HalfEdge = { toVertex, face, edge, opposite, next }
747 // Taking care of new half-edges
748 he0.toVertex = hei_next.toVertex;
749 he0.face = hei.face;
750 he0.edge = new_edgei;
751 he0.opposite = new_hei + 1;
752 he0.next = hei_next.next;
753 he1.toVertex = new_vtxi;
754 he1.face = new_facei;
755 he1.edge = new_edgei;
756 he1.opposite = new_hei;
757 he1.next = new_hei + 2;
758 he2.toVertex = hei.toVertex;
759 he2.face = new_facei;
760 he2.edge = new_edgei + 1;
761 he2.opposite = j;
762 he2.next = hei.next;
763 he3.toVertex = hej_next.toVertex;
764 he3.face = hej.face;
765 he3.edge = new_edgei + 2;
766 he3.opposite = new_hei + 4;
767 he3.next = hej_next.next;
768 he4.toVertex = new_vtxi;
769 he4.face = new_facei + 1;
770 he4.edge = new_edgei + 2;
771 he4.opposite = new_hei + 3;
772 he4.next = new_hei + 5;
773 he5.toVertex = hej.toVertex;
774 he5.face = new_facei + 1;
775 he5.edge = hei.edge;
776 he5.opposite = i;
777 he5.next = hej.next;
778 // Updating existing half-edges
779 hei.toVertex = new_vtxi;
780 hei.opposite = new_hei + 5;
781 hei.next = new_hei;
782 hej.toVertex = new_vtxi;
783 hej.edge = new_edgei + 1;
784 hej.opposite = new_hei + 2;
785 hej.next = new_hei + 3;
786 hei_next.face = new_facei;
787 hei_next.next = new_hei + 1;
788 hej_next.face = new_facei + 1;
789 hej_next.next = new_hei + 4;
790 // Updating other arrays
791 myEdgeHalfEdges.push_back( new_hei );
792 myEdgeHalfEdges.push_back( j );
793 myEdgeHalfEdges.push_back( new_hei + 3 );
794 myFaceHalfEdges.push_back( new_hei + 1 );
795 myFaceHalfEdges.push_back( new_hei + 4 );
796 myFaceHalfEdges[ hei.face ] = i;
797 myFaceHalfEdges[ hej.face ] = j;
798 myEdgeHalfEdges[ hei.edge ] = i;
799 if ( update_arc2index )
800 {
801 const VertexIndex vi = he5.toVertex;
802 const VertexIndex vk = hei_next.toVertex;
803 const VertexIndex vj = he2.toVertex;
804 const VertexIndex vl = hej_next.toVertex;
805 myArc2Index.erase( Arc( vi, vj ) );
806 myArc2Index.erase( Arc( vj, vi ) );
807 myArc2Index[ Arc( vi, new_vtxi ) ] = i;
808 myArc2Index[ Arc( new_vtxi, vi ) ] = new_hei + 5;
809 myArc2Index[ Arc( vj, new_vtxi ) ] = j;
810 myArc2Index[ Arc( new_vtxi, vj ) ] = new_hei + 2;
811 myArc2Index[ Arc( vk, new_vtxi ) ] = new_hei + 1;
812 myArc2Index[ Arc( new_vtxi, vk ) ] = new_hei;
813 myArc2Index[ Arc( vl, new_vtxi ) ] = new_hei + 4;
814 myArc2Index[ Arc( new_vtxi, vl ) ] = new_hei + 3;
815 }
816 return new_vtxi;
817 }
818
827 bool isMergeable( const Index hei ) const {
828 ASSERT( hei != HALF_EDGE_INVALID_INDEX );
829 const HalfEdge& he = halfEdge( hei );
830 const Index hei2 = he.opposite;
831 const HalfEdge& he2 = halfEdge( hei2 );
832 // check if hei borders an infinite face.
833 if ( ( he.face == HALF_EDGE_INVALID_INDEX )
834 || ( he2.face == HALF_EDGE_INVALID_INDEX ) )
835 return false;
836 // Checks that he1 and he2 border a triangle.
837 if ( ( nbSides( hei ) != 3 ) || ( nbSides( hei2 ) != 3) )
838 return false;
839 // Checks that no vertices around lie on the boundary.
840 return ( ! isVertexBoundary( he.toVertex ) )
841 && ( ! isVertexBoundary( he2.toVertex ) )
842 && ( ! isVertexBoundary( halfEdge( he.next ).toVertex ) )
843 && ( ! isVertexBoundary( halfEdge( he2.next ).toVertex ) );
844 }
845
846
865 VertexIndex merge( const Index hei, bool update_arc2index = true ) {
866 const Index i1 = hei; // arc (v1,v2)
867 const HalfEdge& he1 = halfEdge( i1 );
868 const Index i1n = he1.next;
869 const HalfEdge& he1n = halfEdge( i1n );
870 const Index i1nn = he1n.next;
871 const Index i2 = he1.opposite; // arc (v2,v1)
872 const HalfEdge& he2 = halfEdge( i2 );
873 const Index i2n = he2.next;
874 const HalfEdge& he2n = halfEdge( i2n );
875 const Index i2nn = he2n.next;
876 const Index iext1nn = halfEdge( i1nn ).opposite;
877 const Index iext1n = halfEdge( i1n ).opposite;
878 const Index iext2nn = halfEdge( i2nn ).opposite;
879 const Index iext2n = halfEdge( i2n ).opposite;
880 const VertexIndex v1 = he2.toVertex;
881 // Storing v2 the deleted vertex.
882 const VertexIndex v2 = he1.toVertex;
883 const VertexIndex v3 = he1n.toVertex;
884 const VertexIndex v4 = he2n.toVertex;
885 // Storing deleted face indices f1 and f2
886 const FaceIndex f1 = he1.face;
887 const FaceIndex f2 = he2.face;
888 // Storing deleted edge indices (v1,v2), (v2,v3) and (v2,v3').
889 const EdgeIndex e1 = he1.edge;
890 const EdgeIndex e2 = he1n.edge;
891 const EdgeIndex e3 = halfEdge( i2nn ).edge;
892 const EdgeIndex ev13 = halfEdge( iext1nn ).edge;
893 const EdgeIndex ev14 = halfEdge( iext2n ).edge;
894 // For debug
895 auto nbV1 = nbNeighboringVertices( v1 );
896 auto nbV2 = nbNeighboringVertices( v2 );
897 auto nbV3 = nbNeighboringVertices( v3 );
898 auto nbV4 = nbNeighboringVertices( v4 );
899 // Changes toVertex field of half-edges pointing to v2
900 std::vector<VertexIndex> outer_v; // stores vertices around v2
901 std::vector<Index> inner_he; // stores half-edges from v2
902 std::vector<Index> outer_he; // stores half-edges toward v2
903 Index i = i1n;
904 do {
905 const HalfEdge& he = halfEdge( i );
906 const Index iopp = he.opposite;
907 outer_v. push_back( he.toVertex );
908 inner_he.push_back( i );
909 outer_he.push_back( iopp );
910 HalfEdge& heopp = myHalfEdges[ iopp ];
911 ASSERT( heopp.toVertex == v2 );
912 heopp.toVertex = v1;
913 i = heopp.next;
914 } while ( i != i1n ); // i2 precedes i1n around the vertex v2.
915 // std::cout << "#outer_v=" << outer_v.size() << std::endl;
916 // Gluing arcs around the two triangles
917 // std::cout << "Gluing arcs around the two triangles" << std::endl;
918 myHalfEdges[ iext1nn ].opposite = iext1n;
919 myHalfEdges[ iext1n ].opposite = iext1nn;
920 myHalfEdges[ iext2nn ].opposite = iext2n;
921 myHalfEdges[ iext2n ].opposite = iext2nn;
922 // Changing edges of merged edges
923 // std::cout << "Changing edges of merged edges" << std::endl;
924 myHalfEdges[ iext1n ].edge = ev13;
925 myHalfEdges[ iext2nn ].edge = ev14;
926 // Taking care of look-up tables.
927 // (1) myVertexHalfEdges
928 // std::cout << "(1) myVertexHalfEdges" << std::endl;
929 myVertexHalfEdges[ v1 ] = iext1nn;
930 if ( myVertexHalfEdges[ v3 ] == i1nn )
931 myVertexHalfEdges[ v3 ] = iext1n;
932 if ( myVertexHalfEdges[ v4 ] == i2nn )
933 myVertexHalfEdges[ v4 ] = iext2n;
934 // (2) myFaceHalfEdges -> nothing to do.
935 // (3) myEdgeHalfEdges
936 // std::cout << "(3) myEdgeHalfEdges" << std::endl;
937 myEdgeHalfEdges[ ev13 ] = iext1nn;
938 myEdgeHalfEdges[ ev14 ] = iext2n;
939 // (4) myArc2Index only if asked
940 if ( update_arc2index ) {
941 for ( std::size_t j = 0; j < outer_v.size(); ++j ) {
942 myArc2Index.erase( Arc( v2, outer_v[ j ] ) );
943 myArc2Index.erase( Arc( outer_v[ j ], v2 ) );
944 }
945 for ( std::size_t j = 1; j < ( outer_v.size() - 2 ); ++j ) {
946 myArc2Index[ Arc( v1, outer_v[ j ] ) ] = inner_he[ j ];
947 myArc2Index[ Arc( outer_v[ j ], v1 ) ] = outer_he[ j ];
948 }
949 myArc2Index[ Arc( v3, v1 ) ] = iext1n;
950 myArc2Index[ Arc( v1, v4 ) ] = iext2nn;
951 }
952 // Debug
953 auto new_nbV1 = nbNeighboringVertices( v1 );
954 auto new_nbV3 = nbNeighboringVertices( v3 );
955 auto new_nbV4 = nbNeighboringVertices( v4 );
956 if ( ( new_nbV1 != nbV1+nbV2-4 )
957 || ( new_nbV3 != nbV3 -1 )
958 || ( new_nbV4 != nbV4 -1 ) )
959 trace.warning() << "Invalid nb of neighbors: "
960 << " nbV1=" << nbV1
961 << " nbV2=" << nbV2
962 << " nbV3=" << nbV3
963 << " nbV4=" << nbV4 << std::endl
964 << " new_nbV1=" << new_nbV1
965 << " new_nbV3=" << new_nbV3
966 << " new_nbV4=" << new_nbV4 << std::endl;
967
968 // Renumbering of 1 vertex, 3 edges, 2 faces, 6 half-edges
969 renumberVertex( v2, update_arc2index );
970 std::array< EdgeIndex, 3 > E = { e1, e2, e3 };
971 std::sort( E.begin(), E.end(), std::greater<EdgeIndex>() );
972 for ( Index e : E ) {
973 renumberEdge( e );
974 }
975 std::array< FaceIndex, 2 > F = { f1, f2 };
976 std::sort( F.begin(), F.end(), std::greater<FaceIndex>() );
977 for ( Index f : F ) {
978 renumberFace( f );
979 }
980 std::array< Index, 6 > T = { i1, i1n, i1nn, i2, i2n, i2nn };
981 std::sort( T.begin(), T.end(), std::greater<Index>() );
982 for ( Index t : T ) {
983 renumberHalfEdge( t );
984 }
985 return v1;
986 }
987
988 // ------------------------ protected services -------------------------
989 protected:
990
994 void renumberVertex( const VertexIndex vi, bool update_arc2index = true ) {
995 // Vertex j becomes vertex i
996 const VertexIndex vj = nbVertices() - 1;
997 if ( vi != vj ) {
998 const Index j = myVertexHalfEdges[ vj ];
999 Index k = j;
1000 // Turns around vertex vj to modify toVertex fields.
1001 do {
1002 const HalfEdge& hek = halfEdge( k );
1003 const Index kopp = hek.opposite;
1004 HalfEdge& hekopp = myHalfEdges[ kopp ];
1005 ASSERT( hekopp.toVertex == vj );
1006 hekopp.toVertex = vi;
1007 if ( update_arc2index ) {
1008 myArc2Index.erase( Arc( vj, hek.toVertex ) );
1009 myArc2Index.erase( Arc( hek.toVertex, vj ) );
1010 myArc2Index[ Arc( vi, hek.toVertex ) ] = k;
1011 myArc2Index[ Arc( hek.toVertex, vi ) ] = kopp;
1012 }
1013 k = hekopp.next;
1014 } while ( k != j );
1015 myVertexHalfEdges[ vi ] = j;
1016 }
1017 myVertexHalfEdges.pop_back();
1018 }
1019
1023 void renumberEdge( const EdgeIndex ei ) {
1024 const EdgeIndex ej = nbEdges() - 1;
1025 if ( ei != ej ) {
1026 const Index j = myEdgeHalfEdges[ ej ];
1027 HalfEdge& hej = myHalfEdges[ j ];
1028 const Index jopp = hej.opposite;
1029 HalfEdge& hejopp = myHalfEdges[ jopp ];
1030 ASSERT( hej.edge == ej );
1031 ASSERT( hejopp.edge == ej );
1032 hej.edge = ei;
1033 hejopp.edge = ei;
1034 myEdgeHalfEdges[ ei ] = j;
1035 }
1036 myEdgeHalfEdges.pop_back();
1037 }
1038
1042 void renumberFace( const FaceIndex fi ) {
1043 const FaceIndex fj = nbFaces() - 1;
1044 if ( fi != fj ) {
1045 const Index j = myFaceHalfEdges[ fj ];
1046 Index k = j;
1047 do {
1048 HalfEdge& hek = myHalfEdges[ k ];
1049 ASSERT( hek.face == fj );
1050 hek.face = fi;
1051 k = hek.next;
1052 } while ( k != j );
1053 myFaceHalfEdges[ fi ] = j;
1054 }
1055 myFaceHalfEdges.pop_back();
1056 }
1057
1061 void renumberHalfEdge( const Index i, bool update_arc2index = true ) {
1062 const Index j = nbHalfEdges() - 1;
1063 if ( i != j ) {
1064 // std::cout << "renumberHalfEdge j=" << j << std::endl;
1065 const HalfEdge& hej = halfEdge( j );
1066 const Index jopp = hej.opposite;
1067 // std::cout << "renumberHalfEdge jopp=" << jopp << std::endl;
1068 HalfEdge& hejopp = myHalfEdges[ jopp ];
1069 const VertexIndex vj = hejopp.toVertex; // hej corresponds to (vj,vk)
1070 const VertexIndex vk = hej.toVertex; // hejopp corresponds to (vk,vj)
1071 // Update opposite and previous
1072 Index k = hej.next;
1073 while ( halfEdge( k ).next != j ) k = halfEdge( k ).next;
1074 myHalfEdges[ k ].next = i;
1075 hejopp.opposite = i;
1076 // Take care of look-up tables
1077 // std::cout << "myVertexHalfEdges[" << vj << "]" << std::endl;
1078 if ( myVertexHalfEdges[ vj ] == j )
1079 myVertexHalfEdges[ vj ] = i;
1080 // std::cout << "myEdgeHalfEdges[" << hej.edge << "]" << std::endl;
1081 if ( myEdgeHalfEdges[ hej.edge ] == j )
1082 myEdgeHalfEdges[ hej.edge ] = i;
1083 // std::cout << "myFaceHalfEdges[" << hej.face << "]" << std::endl;
1084 if ( myFaceHalfEdges[ hej.face ] == j )
1085 myFaceHalfEdges[ hej.face ] = i;
1086 // std::cout << "myArc2Index[" << vj << "," << vk << "]" << std::endl;
1087 if ( update_arc2index ) {
1088 myArc2Index[ Arc( vj, vk ) ] = i;
1089 }
1090 // Copy last half-edge into i-th half-edge.
1091 myHalfEdges[ i ] = hej;
1092 }
1093 myHalfEdges.pop_back();
1094 }
1095
1096 // ------------------------- consistency services --------------------
1097 public:
1109 bool isValid( bool check_arc2index = true ) const
1110 {
1111 bool ok = true;
1112 // Checks that indices are within range
1113 for ( Index i = 0; i < nbHalfEdges(); i++ )
1114 {
1115 const HalfEdge& he = myHalfEdges[ i ];
1116 if ( he.next >= nbHalfEdges() ) {
1117 trace.warning() << "[HalfEdgeDataStructure::isValid] "
1118 << "half-edge " << i
1119 << " has invalid next half-edge " << he.next
1120 << std::endl;
1121 ok = false;
1122 }
1123 if ( he.opposite >= nbHalfEdges() ) {
1124 trace.warning() << "[HalfEdgeDataStructure::isValid] "
1125 << "half-edge " << i
1126 << " has invalid opposite half-edge " << he.opposite
1127 << std::endl;
1128 ok = false;
1129 }
1130 if ( he.toVertex >= nbVertices() ) {
1131 trace.warning() << "[HalfEdgeDataStructure::isValid] "
1132 << "half-edge " << i
1133 << " has invalid toVertex " << he.toVertex
1134 << std::endl;
1135 ok = false;
1136 }
1137 if ( he.face >= nbFaces() && he.face != HALF_EDGE_INVALID_INDEX ) {
1138 trace.warning() << "[HalfEdgeDataStructure::isValid] "
1139 << "half-edge " << i
1140 << " has invalid face " << he.face
1141 << std::endl;
1142 ok = false;
1143 }
1144 }
1145 // Checks that opposite are correct
1146 for ( Index i = 0; i < nbHalfEdges(); i++ )
1147 {
1148 const Index j = myHalfEdges[ i ].opposite;
1149 if ( j == HALF_EDGE_INVALID_INDEX ) {
1150 trace.warning() << "[HalfEdgeDataStructure::isValid] "
1151 << "half-edge " << i << " has invalid opposite half-edge." << std::endl;
1152 ok = false;
1153 }
1154 if ( myHalfEdges[ j ].opposite != i ) {
1155 trace.warning() << "[HalfEdgeDataStructure::isValid] "
1156 << "half-edge " << i << " has opposite half-edge " << j
1157 << " but the latter has opposite half-edge " << myHalfEdges[ j ].opposite << std::endl;
1158 ok = false;
1159 }
1160 const VertexIndex vi = myHalfEdges[ i ].toVertex;
1161 const VertexIndex vj = myHalfEdges[ j ].toVertex;
1162 if ( vi == vj ) {
1163 trace.warning() << "[HalfEdgeDataStructure::isValid] "
1164 << "half-edge " << i << " and its opposite half-edge " << j
1165 << " have the same toVertex " << vi << std::endl;
1166 ok = false;
1167 }
1168 if ( vi >= nbVertices() ) {
1169 trace.warning() << "[HalfEdgeDataStructure::isValid] "
1170 << "half-edge " << i
1171 << " points to an invalid vertex " << vi << std::endl;
1172 ok = false;
1173 }
1174 if ( vj >= nbVertices() ) {
1175 trace.warning() << "[HalfEdgeDataStructure::isValid] "
1176 << "opposite half-edge " << j
1177 << " points to an invalid vertex " << vj << std::endl;
1178 ok = false;
1179 }
1180 }
1181 // Checks that vertices have a correct starting half-edge.
1182 for ( VertexIndex i = 0; i < nbVertices(); i++ )
1183 {
1184 const Index j = myVertexHalfEdges[ i ];
1185 const Index jopp = myHalfEdges[ j ].opposite;
1186 const VertexIndex ip = myHalfEdges[ jopp ].toVertex;
1187 if ( ip != i ) {
1188 trace.warning() << "[HalfEdgeDataStructure::isValid] "
1189 << "vertex " << i << " is associated to half-edge " << j
1190 << " but its opposite half-edge " << jopp << " points to vertex " << ip << std::endl;
1191 ok = false;
1192 }
1193 }
1194 // Checks that faces have a correct bordering half-edge.
1195 for ( FaceIndex f = 0; f < nbFaces(); f++ )
1196 {
1197 const Index i = myFaceHalfEdges[ f ];
1198 if ( myHalfEdges[ i ].face != f ) {
1199 trace.warning() << "[HalfEdgeDataStructure::isValid] "
1200 << "face " << f << " is associated to half-edge " << i
1201 << " but its associated face is " << myHalfEdges[ i ].face << std::endl;
1202 ok = false;
1203 }
1204 }
1205 // Checks that following next half-edges turns around the same face.
1206 for ( FaceIndex f = 0; f < nbFaces(); f++ )
1207 {
1208 Index i = myFaceHalfEdges[ f ];
1209 const Index start = i;
1210 do {
1211 const HalfEdge& he = halfEdge( i );
1212 if ( he.face != f ) {
1213 trace.warning() << "[HalfEdgeDataStructure::isValid] "
1214 << "when turning around face " << f << ", half-edge " << i
1215 << " is associated to face " << he.face << std::endl;
1216 ok = false;
1217 }
1218 i = he.next;
1219 }
1220 while ( i != start );
1221 }
1222 // Checks that turning around a vertex gives half-edges associated to this vertex.
1223 for ( VertexIndex v = 0; v < nbVertices(); v++ )
1224 {
1225 const Index s = halfEdgeIndexFromVertexIndex( v );
1226 Index i = s;
1227 do
1228 {
1229 const HalfEdge& he = halfEdge( i );
1230 const VertexIndex w = halfEdge( he.opposite ).toVertex;
1231 if ( v != w ) {
1232 trace.warning() << "[HalfEdgeDataStructure::isValid] "
1233 << "when turning around vertex " << v << ", some opposite half-edge "
1234 << he.opposite << " points to " << w << std::endl;
1235 ok = false;
1236 }
1237 i = halfEdge( he.opposite ).next;
1238 }
1239 while ( i != s );
1240 }
1241
1242 // Checks that boundary vertices have specific associated half-edges.
1244 for ( VertexIndex i : bdryV ) {
1245 const Index j = myVertexHalfEdges[ i ];
1246 if ( halfEdge( j ).face != HALF_EDGE_INVALID_INDEX ) {
1247 trace.warning() << "[HalfEdgeDataStructure::isValid] "
1248 << "boundary vertex " << i << " is associated to the half-edge " << j
1249 << " that does not lie on the boundary but on face " << halfEdge( j ).face
1250 << std::endl;
1251 ok = false;
1252 }
1253 }
1254
1255 // Checks that that we can find arcs.
1256 for ( Index i = 0; i < nbHalfEdges(); i++ ) {
1257 const HalfEdge& he = halfEdge( i );
1258 const VertexIndex v2 = he.toVertex;
1259 const VertexIndex v1 = halfEdge( he.opposite ).toVertex;
1260 const Index j = findHalfEdgeIndexFromArc( v1, v2 );
1261 const HalfEdge& he2 = halfEdge( j );
1262 if ( he2.toVertex != v2 ) {
1263 trace.warning() << "[HalfEdgeDataStructure::isValid] "
1264 << "arc (" << v1 << "," << v2 << ") was found to be half-edge " << j
1265 << " but it should be half-edge " << i << std::endl;
1266 ok = false;
1267 }
1268 }
1269
1270 // Checks that edges have two correct associated half-edges.
1271 for ( EdgeIndex ei = 0; ei < nbEdges(); ei++ ) {
1272 const Index i = myEdgeHalfEdges[ ei ];
1273 const HalfEdge& hei = halfEdge( i );
1274 const Index j = hei.opposite;
1275 const HalfEdge& hej = halfEdge( j );
1276 if ( hei.edge != ei ) {
1277 trace.warning() << "[HalfEdgeDataStructure::isValid] "
1278 << "edge " << ei << " is associated to half-edge " << i
1279 << " but its edge is " << hei.edge << std::endl;
1280 ok = false;
1281 }
1282 if ( hej.edge != ei ) {
1283 trace.warning() << "[HalfEdgeDataStructure::isValid] "
1284 << "edge " << ei << " is associated to half-edge " << i
1285 << " of opposite half-edge " << j
1286 << " but its edge is " << hej.edge << std::endl;
1287 ok = false;
1288 }
1289 }
1290
1291 // Checks that arcs have a correct associated half-edge.
1292 if ( check_arc2index )
1293 {
1294 for ( auto arc2idx : myArc2Index )
1295 {
1296 const VertexIndex v1 = arc2idx.first.first;
1297 const VertexIndex v2 = arc2idx.first.second;
1298 const Index i = arc2idx.second;
1299 if ( myHalfEdges[ i ].toVertex != v2 ) {
1300 trace.warning() << "[HalfEdgeDataStructure::isValid] "
1301 << "arc (" << v1 << "," << v2 << ") is associated to half-edge " << i
1302 << " but it points to vertex " << myHalfEdges[ i ].toVertex << std::endl;
1303 ok = false;
1304 }
1305 const Index i2 = myHalfEdges[ i ].opposite;
1306 if ( myHalfEdges[ i2 ].toVertex != v1 ) {
1307 trace.warning() << "[HalfEdgeDataStructure::isValid] "
1308 << "arc (" << v1 << "," << v2 << ") is associated to half-edge " << i
1309 << " but its opposite half-edge " << i2 << " points to vertex " << myHalfEdges[ i2 ].toVertex
1310 << std::endl;
1311 ok = false;
1312 }
1313 }
1314 }
1315 return ok;
1316 }
1317
1323 {
1324 bool ok = true;
1325 // Checks that indices are within range
1326 // Checks that following next half-edges turns around the same face.
1327 for ( FaceIndex f = 0; f < nbFaces(); f++ )
1328 {
1329 Index i = myFaceHalfEdges[ f ];
1330 const Index start = i;
1331 int nbv = 0;
1332 std::set<VertexIndex> V;
1333 std::set<EdgeIndex> E;
1334 do {
1335 const HalfEdge& he = halfEdge( i );
1336 if ( he.face != f ) {
1337 trace.warning() << "[HalfEdgeDataStructure::isValidTriangulation] "
1338 << "when turning around face " << f
1339 << ", half-edge " << i
1340 << " is associated to face " << he.face
1341 << std::endl;
1342 ok = false;
1343 }
1344 V.insert( he.toVertex );
1345 E.insert( he.edge );
1346 nbv++;
1347 i = he.next;
1348 }
1349 while ( i != start );
1350 if ( nbv != 3 ) {
1351 trace.warning() << "[HalfEdgeDataStructure::isValidTriangulation] "
1352 << "when turning around face " << f
1353 << ", we had to visit " << nbv
1354 << " half-edges to loop (instead of 3)" << std::endl;
1355 ok = false;
1356 }
1357 if ( V.size() != 3 ) {
1358 trace.warning() << "[HalfEdgeDataStructure::isValidTriangulation] "
1359 << "when turning around face " << f
1360 << ", the set of vertices has not a size 3:";
1361 for ( auto v : V ) trace.warning() << " " << v;
1362 trace.warning() << std::endl;
1363 ok = false;
1364 }
1365 if ( E.size() != 3 ) {
1366 trace.warning() << "[HalfEdgeDataStructure::isValidTriangulation] "
1367 << "when turning around face " << f
1368 << ", the set of edges has not a size 3:";
1369 for ( auto e : E ) trace.warning() << " " << e;
1370 trace.warning() << std::endl;
1371 ok = false;
1372 }
1373 }
1374 return ok;
1375 }
1376 protected:
1377
1379 std::vector< HalfEdge > myHalfEdges;
1383 std::vector< Index > myVertexHalfEdges;
1387 std::vector< Index > myFaceHalfEdges;
1391 std::vector< Index > myEdgeHalfEdges;
1394
1395
1396 // ----------------------- Interface --------------------------------------
1397 public:
1398
1403 void selfDisplay ( std::ostream & out ) const;
1404
1405 // ------------------------- Protected Datas ------------------------------
1406 private:
1407 // ------------------------- Private Datas --------------------------------
1408 private:
1409 // ------------------------- Hidden services ------------------------------
1410 protected:
1411
1412 static
1414 VertexIndex vi, VertexIndex vj )
1415 {
1416 ASSERT( !de2fi.empty() );
1417
1418 Arc2FaceIndex::const_iterator it = de2fi.find( Arc( vi, vj ) );
1419 // If no such directed edge exists, then there's no such face in the mesh.
1420 // The edge must be a boundary edge.
1421 // In this case, the reverse orientation edge must have a face.
1422 if( it == de2fi.end() )
1423 {
1424 ASSERT( de2fi.find( Arc( vj, vi ) ) != de2fi.end() );
1426 }
1427 return it->second;
1428 }
1429
1430 }; // end of class HalfEdgeDataStructure
1431
1432
1439 std::ostream&
1440 operator<< ( std::ostream & out, const HalfEdgeDataStructure & object );
1441
1442} // namespace DGtal
1443
1444
1446// Includes inline functions.
1447#include "DGtal/topology/HalfEdgeDataStructure.ih"
1448
1449// //
1451
1452#endif // !defined HalfEdgeDataStructure_h
1453
1454#undef HalfEdgeDataStructure_RECURSES
1455#endif // else defined(HalfEdgeDataStructure_RECURSES)
Aim: This class represents an half-edge data structure, which is a structure for representing the top...
std::vector< VertexIndex > VertexIndexRange
Index HalfEdgeIndex
The type used for numbering half-edges (alias)
Arc2Index myArc2Index
The mapping between arcs to their half-edge index.
bool build(const std::vector< PolygonalFace > &polygonal_faces)
const HalfEdge & halfEdge(const Index i) const
Index halfEdgeIndexFromFaceIndex(const FaceIndex fi) const
Index halfEdgeIndexFromEdgeIndex(const EdgeIndex ei) const
void clear()
Clears the data structure.
void selfDisplay(std::ostream &out) const
FaceIndexRange neighboringFaces(const VertexIndex vi) const
VertexIndex split(const Index i, bool update_arc2index=true)
Index findHalfEdgeIndexFromArc(const VertexIndex i, const VertexIndex j) const
static Size getUnorderedEdgesFromPolygonalFaces(const std::vector< PolygonalFace > &polygonal_faces, std::vector< Edge > &edges_out)
std::vector< VertexIndex > PolygonalFace
std::map< Arc, FaceIndex > Arc2FaceIndex
bool isMergeable(const Index hei) const
void renumberFace(const FaceIndex fi)
void renumberEdge(const EdgeIndex ei)
static Size getUnorderedEdgesFromTriangles(const std::vector< Triangle > &triangles, std::vector< Edge > &edges_out)
void getNeighboringFaces(const VertexIndex vi, FaceIndexRange &result) const
std::vector< Arc > boundaryArcs() const
Index FaceIndex
The type for numbering faces.
static FaceIndex arc2FaceIndex(const Arc2FaceIndex &de2fi, VertexIndex vi, VertexIndex vj)
void renumberVertex(const VertexIndex vi, bool update_arc2index=true)
bool isValid(bool check_arc2index=true) const
std::vector< Index > myVertexHalfEdges
VertexIndexRange verticesOfFace(const FaceIndex f) const
Arc arcFromHalfEdgeIndex(const Index i) const
Index EdgeIndex
The type for numbering edges.
Index findHalfEdgeIndexFromArc(const Arc &arc) const
std::vector< HalfEdgeIndex > HalfEdgeIndexRange
bool isFlippable(const Index hei) const
std::size_t Index
The type used for numbering half-edges (an offset an the half-edges structure).
VertexIndexRange neighboringVertices(const VertexIndex vi) const
bool isVertexBoundary(const VertexIndex vi) const
Index VertexIndex
The type for numbering vertices.
std::vector< FaceIndex > FaceIndexRange
Size nbSidesOfFace(const FaceIndex f) const
std::vector< HalfEdge > myHalfEdges
Stores all the half-edges.
bool build(const Size num_vertices, const std::vector< PolygonalFace > &polygonal_faces, const std::vector< Edge > &edges)
bool build(const std::vector< Triangle > &triangles)
std::size_t Size
The type for counting elements.
VertexIndexRange boundaryVertices() const
void renumberHalfEdge(const Index i, bool update_arc2index=true)
std::vector< Index > boundaryHalfEdgeIndices() const
HalfEdgeDataStructure()
Default constructor. The data structure is empty.
void flip(const Index hei, bool update_arc2index=true)
std::vector< EdgeIndex > EdgeIndexRange
std::pair< VertexIndex, VertexIndex > Arc
An arc is a directed edge from a first vertex to a second vertex.
Index halfEdgeIndexFromArc(const Arc &arc) const
VertexIndex merge(const Index hei, bool update_arc2index=true)
Index halfEdgeIndexFromArc(const VertexIndex i, const VertexIndex j) const
Index halfEdgeIndexFromVertexIndex(const VertexIndex vi) const
Size nbNeighboringVertices(const Index vi) const
bool build(const Size num_vertices, const std::vector< Triangle > &triangles, const std::vector< Edge > &edges)
std::map< Arc, Index > Arc2Index
void getNeighboringVertices(const VertexIndex vi, VertexIndexRange &result) const
std::ostream & warning()
DGtal is the top-level namespace which contains all DGtal functions and types.
static std::size_t const HALF_EDGE_INVALID_INDEX
std::ostream & operator<<(std::ostream &out, const ClosedIntegerHalfPlane< TSpace > &object)
Trace trace
Definition: Common.h:154
Edge(VertexIndex vi, VertexIndex vj)
bool operator<(const Edge &other) const
Index next
Index into the halfedges array.
HalfEdge()
Default constructor. The half-edge is invalid.
EdgeIndex edge
Index into the edges array.
VertexIndex toVertex
The end vertex of this half-edge as an index into the vertex array.
Index opposite
Index into the halfedges array.
FaceIndex face
Index into the face array.
Represents an unoriented triangle as three vertices.
Triangle(VertexIndex v0, VertexIndex v1, VertexIndex v2)
std::array< VertexIndex, 3 > v
The three vertex indices.