Add missing constraint to flips in HalfEdgeDataStructure and TriangulatedSurface

ChangeLog.md

@@ -128,7 +128,7 @@
     Shortcuts::makePrimalPolygonalSurface, as well as
     MeshHelpers::digitalSurface2PrimalPolygonalSurface, all use the
     CCW ordering by default (in 3D).
-    (Jacques-Olivier Lachaud,[#1421](https://github.com/DGtal-team/DGtal/pull/1445))
+    (Jacques-Olivier Lachaud,[#1445](https://github.com/DGtal-team/DGtal/pull/1445))
 
 - *images*
   - Fix the image origin that was not taken into account in class
@@ -137,18 +137,23 @@
   - Add domainShift to ImageContainerByITKImage.
     (Pablo Hernandez-Cerdan,
     [#1490](https://github.com/DGtal-team/DGtal/pull/1490))
-    
+
 - *IO*
   - Removing a `using namespace std;` in the Viewer3D hearder file. (David
     Coeurjolly [#1413](https://github.com/DGtal-team/DGtal/pull/1413))
   - Fixing cast from const to mutable iterator in GradientColorMap.
     (Roland Denis [#1486](https://github.com/DGtal-team/DGtal/pull/1486))
 
+- *Topology*
+  - Add missing constraint to flips in HalfEdgeDataStructure
+    (Jacques-Olivier Lachaud,[#1498](https://github.com/DGtal-team/DGtal/pull/1498))
 
-- *Shapes package*
+- *Shapes*
   - Fix bug in Astroid parameter() method : orientation correction
    (Adrien Krähenbühl,
    [#1325](https://github.com/DGtal-team/DGtal/pull/1426))
+  - Add missing constraint to flips in TriangulatedSurface
+    (Jacques-Olivier Lachaud,[#1498](https://github.com/DGtal-team/DGtal/pull/1498))
 
 - *DEC*
   - Fix issue (https://github.com/DGtal-team/DGtal/issues/1441)

src/DGtal/shapes/TriangulatedSurface.h

@@ -569,8 +569,9 @@ namespace DGtal
   public:
 
     /// An edge is (topologically) flippable iff: (1) it does not lie
-    /// on the boundary, (2) it is bordered by two triangles. (2) is
-    /// always true for a triangulated surface.
+    /// on the boundary, (2) it is bordered by two triangles, (3) the
+    /// two other vertices of the quad are not already neighbors. (2)
+    /// is always true for a triangulated surface.
     ///
     /// @param a any arc.
     /// @return 'true' if the edge containing \a a is topologically flippable.

src/DGtal/shapes/TriangulatedSurface.ih

@@ -441,8 +441,14 @@ inline
 bool
 DGtal::TriangulatedSurface<TPoint>::isFlippable( const Arc a ) const
 {
-  return ( ! isArcBoundary( a ) )
-    && ( ! isArcBoundary( opposite( a ) ) );
+  if ( isArcBoundary( a ) || isArcBoundary( opposite( a ) ) ) return false;
+  const auto v1 = head( next( a ) );
+  const auto v2 = head( next( opposite( a ) ) );
+  std::vector< Index > neighbors_v1;
+  auto outIt = std::back_inserter( neighbors_v1 );
+  writeNeighbors( outIt, v1 );
+  const auto itv2 = std::find( neighbors_v1.cbegin(), neighbors_v1.cend(), v2 );
+  return itv2 == neighbors_v1.cend();
 }
 
 //-----------------------------------------------------------------------------

src/DGtal/topology/HalfEdgeDataStructure.h

@@ -622,7 +622,8 @@ namespace DGtal
   public:
 
     /// An edge is (topologically) flippable iff: (1) it does not lie
-    /// on the boundary, (2) it is bordered by two triangles.
+    /// on the boundary, (2) it is bordered by two triangles, (3) the
+    /// two other vertices of the quad are not already neighbors.
     ///
     /// @param hei any valid half-edge index.
     /// @return 'true' if the edge containing \a hei is topologically flippable.
@@ -631,14 +632,22 @@ namespace DGtal
     bool isFlippable( const Index hei ) const
     {
       ASSERT( hei != HALF_EDGE_INVALID_INDEX );
-      const HalfEdge& he = halfEdge( hei );
-      const Index    hei2 = he.opposite;
+      const HalfEdge&   he = halfEdge( hei );
+      const Index     hei2 = he.opposite;
+      const HalfEdge&  he2 = halfEdge( hei2 );
       // check if hei borders an infinite face.
-      if ( he.face == HALF_EDGE_INVALID_INDEX ) return false; 
+      if ( he.face  == HALF_EDGE_INVALID_INDEX ) return false; 
       // check if hei2 borders an infinite face.
-      if ( halfEdge( hei2 ).face == HALF_EDGE_INVALID_INDEX ) return false; 
+      if ( he2.face == HALF_EDGE_INVALID_INDEX ) return false; 
       // Checks that he1 and he2 border a triangle.
-      return ( nbSides( hei ) == 3 ) && ( nbSides( hei2 ) == 3);
+      if ( ( nbSides( hei ) != 3 ) || ( nbSides( hei2 ) != 3 ) ) return false;
+      // Checks that the two other vertices of the two surrounding
+      // triangles are not already neighbors
+      const VertexIndex v1 = halfEdge( he.next  ).toVertex;
+      const VertexIndex v2 = halfEdge( he2.next ).toVertex;
+      const auto neighb_v1 = neighboringVertices( v1 );
+      const auto     it_v2 = std::find( neighb_v1.cbegin(), neighb_v1.cend(), v2 );
+      return it_v2 == neighb_v1.cend();
     }
 
     /// Tries to flip the edge containing \a hei.

tests/topology/testHalfEdgeDataStructure.cpp

@@ -495,14 +495,14 @@ SCENARIO( "HalfEdgeDataStructure flips", "[halfedge][flips]" ){
   }
   GIVEN( "A tetrahedron" ) {
     HalfEdgeDataStructure mesh = makeTetrahedron();
-    THEN( "All edges are flippable" ) {
+    THEN( "No edges are flippable" ) {
       int nbflippable = 0;
       for ( Size e = 0; e < mesh.nbEdges(); e++ )
         {
           if ( mesh.isFlippable( mesh.halfEdgeIndexFromEdgeIndex( e ) ) )
             nbflippable++;
         }
-      REQUIRE( nbflippable == mesh.nbEdges() );
+      REQUIRE( nbflippable == 0 );
     }
   }
   GIVEN( "Two triangles incident by an edge" ) {