Complete testSurfaceMesh

src/DGtal/shapes/SurfaceMesh.h

@@ -880,6 +880,11 @@ namespace DGtal
  /// @note Time complexity is O(1).
  bool isFlippable( const Edge e ) const;
 
+ /// @pre `isFlippable( e )` must be true.
+ /// @param e any edge.
+ /// @return the two other vertices of the quadrilateral around the edge \a e. 
+ VertexPair otherDiagonal( const Edge e ) const;
+
  /// Flip the edge \a e. Be careful that after the flip, this edge
  /// index determines another edge, which is the other diagonal of
  /// the quadrilateral having \a e as its diagonal.

src/DGtal/shapes/SurfaceMesh.ih

@@ -894,6 +894,29 @@ isFlippable( const Edge e ) const
  return itl == Nk.cend();
 }
 
+//-----------------------------------------------------------------------------
+template <typename TRealPoint, typename TRealVector>
+typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::VertexPair
+DGtal::SurfaceMesh<TRealPoint, TRealVector>::
+otherDiagonal( const Edge e ) const
+{
+ // only valid if `isFlippable( e )` is true.
+ const Faces& rfaces = edgeRightFaces( e );
+ const Faces& lfaces = edgeLeftFaces ( e );
+ const Face rf = rfaces.front(); //< some `(..., j, i, ... )` since faces are ccw.
+ const Face lf = lfaces.front(); //< some `(..., i, j, ... )` since faces are ccw.
+ const Vertices& rvtx = incidentVertices( rf );
+ const Vertices& lvtx = incidentVertices( lf );
+ Vertex i, j;
+ std::tie( i, j ) = edgeVertices( e );
+ const auto ir = ( rvtx[ 0 ] == i ) ? 0 : ( ( rvtx[ 1 ] == i ) ? 1 : 2 );
+ const auto il = ( lvtx[ 0 ] == i ) ? 0 : ( ( lvtx[ 1 ] == i ) ? 1 : 2 );
+ const Vertex k = rvtx[ (ir + 1) % 3 ]; // right triangle is (j,i,k)
+ const Vertex l = lvtx[ (il + 2) % 3 ]; // left triangle is (i,j,l).
+ return ( k < l ) ? std::make_pair( k, l ) : std::make_pair( l, k );
+}
+
+
 //-----------------------------------------------------------------------------
 template <typename TRealPoint, typename TRealVector>
 void
@@ -917,8 +940,8 @@ flip( const Edge e, bool recompute_face_normals )
  // (1) We must collect all information: right and left face, vertices k and l
  const Face rf = edgeRightFaces( e ).front(); //< some `(..., j, i, ... )` since faces are ccw.
  const Face lf = edgeLeftFaces ( e ).front(); //< some `(..., i, j, ... )` since faces are ccw.
- Vertices& rvtx = incidentVertices( rf );
- Vertices& lvtx = incidentVertices( lf );
+ Vertices& rvtx = myIncidentVertices[ rf ];
+ Vertices& lvtx = myIncidentVertices[ lf ];
  Vertex i, j;
  std::tie( i, j ) = edgeVertices( e );
  const auto ir = ( rvtx[ 0 ] == i ) ? 0 : ( ( rvtx[ 1 ] == i ) ? 1 : 2 );

tests/shapes/testSurfaceMesh.cpp

@@ -104,6 +104,27 @@ PointVector<3,double> > makeNonManifoldBoundary()
  faces.cbegin(), faces.cend() );
 }
 
+SurfaceMesh< PointVector<3,double>,
+PointVector<3,double> > makeTetrahedron()
+{
+ typedef PointVector<3,double> RealPoint;
+ typedef PointVector<3,double> RealVector;
+ typedef SurfaceMesh< RealPoint, RealVector > PolygonMesh;
+ typedef PolygonMesh::Vertices Vertices;
+ std::vector< RealPoint > positions;
+ std::vector< Vertices > faces;
+ positions.push_back( RealPoint( 0, 0, 0 ) );
+ positions.push_back( RealPoint( 1, 0, 0 ) );
+ positions.push_back( RealPoint( 0, 1, 0 ) );
+ positions.push_back( RealPoint( 0, 0, 1 ) );
+ faces.push_back( { 0, 1, 2 } );
+ faces.push_back( { 1, 0, 3 } );
+ faces.push_back( { 0, 2, 3 } ); 
+ faces.push_back( { 3, 2, 1 } );
+ return PolygonMesh( positions.cbegin(), positions.cend(),
+ faces.cbegin(), faces.cend() );
+}
+
 SCENARIO( "SurfaceMesh< RealPoint3 > concept check tests", "[surfmesh][concepts]" )
 {
  typedef PointVector<3,double> RealPoint;
@@ -335,9 +356,9 @@ SCENARIO( "SurfaceMesh< RealPoint3 > reader/writer tests", "[surfmesh][io]" )
 
 SCENARIO( "SurfaceMesh< RealPoint3 > boundary tests", "[surfmesh][boundary]" )
 {
- auto polymesh = makeNonManifoldBoundary();
+ auto polymesh  = makeNonManifoldBoundary();
  auto polymesh2 = makeBox();
- WHEN( "Checking the topolopgy of the mesh boundary" ) {
+ WHEN( "Checking the topology of the mesh boundary" ) {
  auto chains = polymesh2.computeManifoldBoundaryChains();
  THEN( "The box as a manifold boundary" ) {
  CAPTURE(chains);
@@ -351,3 +372,143 @@ SCENARIO( "SurfaceMesh< RealPoint3 > boundary tests", "[surfmesh][boundary]" )
  }
  }
 }
+
+SCENARIO( "SurfaceMesh< RealPoint3 > flippable tests", "[surfmesh][flip]" )
+{
+ typedef PointVector<3,double> RealPoint;
+ typedef PointVector<3,double> RealVector;
+ typedef SurfaceMesh< RealPoint, RealVector > PolygonMesh;
+ typedef SurfaceMeshHelper< RealPoint, RealVector > PolygonMeshHelper;
+ typedef PolygonMeshHelper::NormalsType NormalsType;
+ auto meshBox = makeBox();
+ auto meshTetra = makeTetrahedron();
+ auto meshLantern = PolygonMeshHelper::makeLantern( 3.0, 3.0, RealPoint::zero,
+ 10, 10, NormalsType::NO_NORMALS );
+ auto meshTorus = PolygonMeshHelper::makeTorus( 3.0, 1.0, RealPoint::zero,
+ 10, 10, 0, NormalsType::NO_NORMALS );
+ WHEN( "Checking if one can flip box edges" ) {
+ auto nb_flippable = 0;
+ for ( auto e = 0; e < meshBox.nbEdges(); e++ )
+ if ( meshBox.isFlippable( e ) ) nb_flippable++;
+ THEN( "No box edges are flippable (they border quads)" ) {
+ REQUIRE( nb_flippable == 0 );
+ }
+ }
+ WHEN( "Checking if one can flip tetrahedron edges" ) {
+ auto nb_flippable = 0;
+ for ( auto e = 0; e < meshTetra.nbEdges(); e++ )
+ if ( meshTetra.isFlippable( e ) ) nb_flippable++;
+ THEN( "No tetrahedron edges are flippable (the neihgborhood is not simply connected)" ) {
+ REQUIRE( nb_flippable == 0 );
+ }
+ }
+ WHEN( "Checking if one can flip torus edges" ) {
+ auto nb_flippable = 0;
+ for ( auto e = 0; e < meshTorus.nbEdges(); e++ )
+ if ( meshTorus.isFlippable( e ) ) nb_flippable++;
+ THEN( "All torus edges are flippable (it is a closed triangulated surface)" ) {
+ REQUIRE( nb_flippable == meshTorus.nbEdges() );
+ }
+ }
+ WHEN( "Checking if one can flip lantern edges" ) {
+ auto bdry_edges = meshLantern.computeManifoldBoundaryEdges();
+ auto inner_edges = meshLantern.computeManifoldInnerEdges();
+ auto nb_flippable = 0;
+ auto nb_bdry_flippable = 0;
+ auto nb_inner_flippable = 0; 
+ for ( auto e = 0; e < meshLantern.nbEdges(); e++ )
+ if ( meshLantern.isFlippable( e ) ) nb_flippable++;
+ for ( auto e : bdry_edges )
+ if ( meshLantern.isFlippable( e ) ) nb_bdry_flippable++;
+ for ( auto e : inner_edges )
+ if ( meshLantern.isFlippable( e ) ) nb_inner_flippable++;
+ THEN( "Innner lantern edges are flippable while boundary edges are not flippable" ) {
+ REQUIRE( nb_flippable == inner_edges.size() );
+ REQUIRE( nb_bdry_flippable == 0 );
+ REQUIRE( nb_flippable == nb_inner_flippable );
+ REQUIRE( nb_flippable == ( meshLantern.nbEdges() - bdry_edges.size() ) );
+ }
+ }
+}
+
+SCENARIO( "SurfaceMesh< RealPoint3 > flip tests", "[surfmesh][flip]" )
+{
+ typedef PointVector<3,double> RealPoint;
+ typedef PointVector<3,double> RealVector;
+ typedef SurfaceMesh< RealPoint, RealVector > PolygonMesh;
+ typedef SurfaceMeshHelper< RealPoint, RealVector > PolygonMeshHelper;
+ typedef SurfaceMeshWriter< RealPoint, RealVector > PolygonMeshWriter;
+ typedef PolygonMeshHelper::NormalsType NormalsType;
+ auto meshLantern = PolygonMeshHelper::makeLantern( 3.0, 3.0, RealPoint::zero,
+ 10, 10, NormalsType::NO_NORMALS );
+ auto bdry_edges = meshLantern.computeManifoldBoundaryEdges();
+ auto euler = meshLantern.Euler();
+ auto nb_flipped = 0;
+ for ( auto i = 0; i < 100; i++ )
+ {
+ auto e = rand() % meshLantern.nbEdges();
+ if ( meshLantern.isFlippable( e ) )
+ {
+ meshLantern.flip( e, false );
+ nb_flipped++;
+ }
+ }
+ WHEN( "Flipping 100 random edges" ) {
+ THEN( "More than 50 edges were flipped" ) {
+ REQUIRE( nb_flipped > 50 );
+ }
+ THEN( "Euler number is not changed" ) {
+ auto post_euler = meshLantern.Euler();
+ REQUIRE( euler == post_euler );
+ }
+ THEN( "Boundary is unchanged" ) {
+ auto post_bdry_edges = meshLantern.computeManifoldBoundaryEdges(); 
+ REQUIRE( bdry_edges.size() == post_bdry_edges.size() );
+ }
+ }
+}
+
+SCENARIO( "SurfaceMesh< RealPoint3 > restore lantern with flips tests", "[surfmesh][flip]" )
+{
+ typedef PointVector<3,double> RealPoint;
+ typedef PointVector<3,double> RealVector;
+ typedef SurfaceMesh< RealPoint, RealVector > PolygonMesh;
+ typedef SurfaceMeshHelper< RealPoint, RealVector > PolygonMeshHelper;
+ typedef SurfaceMeshWriter< RealPoint, RealVector > PolygonMeshWriter;
+ typedef PolygonMeshHelper::NormalsType NormalsType;
+ auto meshLantern = PolygonMeshHelper::makeLantern( 3.0, 3.0, RealPoint::zero,
+ 10, 10, NormalsType::NO_NORMALS );
+ {
+ std::ofstream output( "lantern.obj" );
+ bool okw = PolygonMeshWriter::writeOBJ( output, meshLantern );
+ output.close();
+ }
+ auto nb_flipped = 0;
+ const auto& X = meshLantern.positions();
+ for ( auto e = 0; e < meshLantern.nbEdges(); e++ )
+ {
+ if ( meshLantern.isFlippable( e ) )
+ {
+ auto ij = meshLantern.edgeVertices ( e );
+ auto kl = meshLantern.otherDiagonal( e );
+ double l2_ij = ( X[ ij.first ] - X[ ij.second ] ).squaredNorm();
+ double l2_kl = ( X[ kl.first ] - X[ kl.second ] ).squaredNorm();
+ if ( l2_kl < l2_ij )
+ {
+ meshLantern.flip( e, false );
+ nb_flipped++;
+ }
+ }
+ }
+ {
+ std::ofstream output( "flipped-lantern.obj" );
+ bool okw = PolygonMeshWriter::writeOBJ( output, meshLantern );
+ output.close();
+ }
+ WHEN( "Flipping all long edges" ) {
+ THEN( "80 edges were flipped" ) {
+ REQUIRE( nb_flipped == 80 );
+ }
+ }
+}
+