SurfaceMesh.h

 
#pragma once
 
#if defined(SurfaceMesh_RECURSES)
#error Recursive header files inclusion detected in SurfaceMesh.h
#else // defined(SurfaceMesh_RECURSES)
#define SurfaceMesh_RECURSES
 
#if !defined SurfaceMesh_h
#define SurfaceMesh_h
 
// Inclusions
#include <iostream>
#include <sstream>
#include <string>
#include "DGtal/base/Common.h"
#include "DGtal/base/IntegerSequenceIterator.h"
#include "DGtal/helpers/StdDefs.h"
 
namespace DGtal
{
  // template class SurfaceMesh
  template < typename TRealPoint, typename TRealVector >
  struct SurfaceMesh
  {
    typedef TRealPoint                              RealPoint;
    typedef TRealVector                             RealVector;
    typedef SurfaceMesh< RealPoint, RealVector >    Self;
    
    static const Dimension dimension = RealPoint::dimension;
    BOOST_STATIC_ASSERT( ( dimension == 3 ) );
    
    typedef typename RealVector::Component          Scalar;
    typedef std::vector<Scalar>                     Scalars;
    typedef std::size_t                             Size;
    typedef std::size_t                             Index;
    typedef Index                                   Face;
    typedef Index                                   Edge;
    typedef Index                                   Vertex;
    typedef std::pair< Edge, Scalar >               WeightedEdge;
    typedef std::pair< Face, Scalar >               WeightedFace;
    typedef std::vector< Vertex >                   Vertices;
    typedef std::vector< Edge >                     Edges;
    typedef std::vector< WeightedEdge >             WeightedEdges;
    typedef std::vector< Face >                     Faces;
    typedef std::vector< WeightedFace >             WeightedFaces;
    typedef std::pair< Vertex, Vertex >             VertexPair;
 
    // Required by CUndirectedSimpleLocalGraph
    typedef std::set<Vertex>                   VertexSet;
    template <typename Value> struct           VertexMap {
      typedef typename std::map<Vertex, Value> Type;
    };
 
    // Required by CUndirectedSimpleGraph
    
    typedef IntegerSequenceIterator< Vertex >       ConstIterator;
    
    //---------------------------------------------------------------------------
  public:
    
    ~SurfaceMesh() = default;
    SurfaceMesh() = default;
    SurfaceMesh( const Self& other ) = default;
    SurfaceMesh( Self&& other ) = default;
    Self& operator=( const Self& other ) = default;
 
    template <typename RealPointIterator, typename VerticesIterator>
    SurfaceMesh( RealPointIterator itPos, RealPointIterator itPosEnd,
                 VerticesIterator itVertices, VerticesIterator itVerticesEnd );
 
    template <typename RealPointIterator, typename VerticesIterator>
    bool init( RealPointIterator itPos, RealPointIterator itPosEnd,
               VerticesIterator itVertices, VerticesIterator itVerticesEnd );
 
    void clear();
 
    
    //---------------------------------------------------------------------------
  public:
    
    template <typename RealVectorIterator>
    bool setVertexNormals( RealVectorIterator itN, RealVectorIterator itNEnd );
 
    template <typename RealVectorIterator>
    bool setFaceNormals( RealVectorIterator itN, RealVectorIterator itNEnd );
 
    void computeFaceNormalsFromPositions();
 
    void computeFaceNormalFromPositions( const Face f );
 
    void computeFaceNormalsFromVertexNormals();
 
    void computeVertexNormalsFromFaceNormals();
    
    void computeVertexNormalsFromFaceNormalsWithMaxWeights();
 
    template <typename AnyRing>
    std::vector<AnyRing> computeFaceValuesFromVertexValues
    ( const std::vector<AnyRing>& vvalues ) const;
    
    template <typename AnyRing>
    std::vector<AnyRing> computeVertexValuesFromFaceValues
    ( const std::vector<AnyRing>& fvalues ) const;
 
    std::vector<RealVector> computeFaceUnitVectorsFromVertexUnitVectors
    ( const std::vector<RealVector>& vuvectors ) const;
 
    std::vector<RealVector> computeVertexUnitVectorsFromFaceUnitVectors
    ( const std::vector<RealVector>& fuvectors ) const;
    
 
    //---------------------------------------------------------------------------
  public:
 
    Size nbVertices() const
    { return myIncidentFaces.size(); }
 
    Size nbEdges() const
    { return myEdgeVertices.size(); }
 
    Size nbFaces() const
    { return myIncidentVertices.size(); }
 
    long Euler() const
    { return nbVertices() - nbEdges() + nbFaces(); }
    
    Edge makeEdge( Vertex i, Vertex j ) const;
 
    const Vertices&  incidentVertices( Face f ) const
    { return myIncidentVertices[ f ]; }
 
    const Faces& incidentFaces( Vertex v ) const
    { return myIncidentFaces[ v ]; }
    
    const Faces& neighborFaces( Face f ) const
    { return myNeighborFaces[ f ]; }
 
    const Vertices& neighborVertices( Vertex v ) const
    { return myNeighborVertices[ v ]; }
 
    VertexPair edgeVertices( Edge e ) const
    { return myEdgeVertices[ e ]; }
    
    const Faces& edgeFaces( Edge e ) const
    { return myEdgeFaces[ e ]; }
 
    const Faces& edgeRightFaces( Edge e ) const 
    { return myEdgeRightFaces[ e ]; }
 
    const Faces& edgeLeftFaces( Edge e ) const 
    { return myEdgeLeftFaces[ e ]; }
 
    const std::vector< Vertices >& allIncidentVertices() const
    { return myIncidentVertices; }
 
    const std::vector< Faces >& allIncidentFaces() const
    { return myIncidentFaces; }
    
    const std::vector< Faces >& allNeighborFaces() const
    { return myNeighborFaces; }
 
    const std::vector< Vertices >& allNeighborVertices() const
    { return myNeighborVertices; }
 
    const std::vector< VertexPair >& allEdgeVertices() const
    { return myEdgeVertices; }
    
    const std::vector< Faces >& allEdgeFaces() const
    { return myEdgeFaces; }
 
    const std::vector< Faces >& allEdgeRightFaces() const 
    { return myEdgeRightFaces; }
 
    const std::vector< Faces >& allEdgeLeftFaces() const 
    { return myEdgeLeftFaces; }
    
 
    //---------------------------------------------------------------------------
  public:
 
    Edges computeManifoldBoundaryEdges() const;
    Edges computeManifoldInnerEdges() const;
    Edges computeManifoldInnerConsistentEdges() const;
    Edges computeManifoldInnerUnconsistentEdges() const;
    Edges computeNonManifoldEdges() const;
    
    bool isBoundariesManifold( bool checkClosed = true ) const
    {
      // computes unordered list of boundary vertices
      std::map<Vertex,Vertices> adjacent;
      auto MBE = this->computeManifoldBoundaryEdges();
      if ( MBE.size() == 0 ) return false;
      
      for (auto e : MBE)
      {
        auto ij = this->edgeVertices(e);
        adjacent[ ij.first ].push_back( ij.second );
        if ( adjacent[ ij.first ] .size() > 2 ) return false;
        adjacent[ ij.second ].push_back( ij.first );
        if ( adjacent[ ij.second ].size() > 2 )  return false;
      }
      //we may check if all curves are closed.
      if ( checkClosed )
        for ( const auto &adj : adjacent )
          if ( adj.second.size() != 2 ) return false;
      return true;
    }
    
    std::vector<Vertices> computeManifoldBoundaryChains() const
    {
      std::vector<Vertices> boundaries;
      Vertices boundary;
      auto MBE = this->computeManifoldBoundaryEdges();
      std::map<Vertex,bool> visited;
      std::map<Vertex,Vertices> adjacent;
      
      ASSERT_MSG(MBE.size()>0,"The surface mesh must have boundary edges");
      ASSERT_MSG(this->isBoundariesManifold(), "The surface mesh mush have manifold boundaries");
      
      //Coompute adjecency relationships
      for (auto e : MBE)
      {
        auto ij = this->edgeVertices(e);
        visited[ij.first] = false;
        visited[ij.second] = false;
        adjacent[ij.first].push_back(ij.second);
        adjacent[ij.second].push_back(ij.first);
      }
      
      auto boundary_it = visited.begin();
      
      while(boundary_it != visited.end() )
      {
        Vertex first = (*boundary_it).first;
        visited[first] = true;
        boundary.clear();
        boundary.push_back(first);
        
        Vertex current = first;
        size_t nb_iter = 0;
        bool pushed=false;
        
        while ((!pushed) && (nb_iter < MBE.size()*2))
        {
          bool ok = false;
          for (auto other : adjacent[current])
            if (!visited[other])
            {
              boundary.push_back(other);
              current = other;
              visited[other] = true;
              ok = true;
              break;
            }
          if (!ok)
          {
            //all neighboors are visited
            for (auto other : adjacent[current])
              if (other == first)
              {
                boundaries.push_back(boundary);
                pushed = true;
                break;
              }
            //if first vertex isn't found then this chain is not
            //homeomorphic to a circle, hence isn't added to boundaries
          }
          nb_iter++;
        }
        boundary_it = std::find_if(visited.begin(), visited.end(),
                                   []
                                   (std::pair<Vertex,bool> x){return !x.second;});
      }
      return boundaries;
    }
    
    // ----------------------- Undirected simple graph services ----------------------
  public:
    
    Size size() const
    { return nbVertices(); }
 
    Size bestCapacity() const
    { return 8; }
      
    Size degree( const Vertex & v ) const
    { return myNeighborVertices[ v ].size(); }
 
    template <typename OutputIterator>
    void
    writeNeighbors( OutputIterator &it ,
                    const Vertex & v ) const
    {
      for ( auto&& nv : myNeighborVertices[ v ] )
        *it++ = nv;
    }
 
    template <typename OutputIterator, typename VertexPredicate>
    void
    writeNeighbors( OutputIterator &it ,
                    const Vertex & v,
                    const VertexPredicate & pred) const
    {
      for ( auto&& nv : myNeighborVertices[ v ] )
        if ( pred( nv ) ) *it++ = nv;
    }
 
    ConstIterator begin() const
    { return ConstIterator( 0 ); }
 
    ConstIterator end() const
    { return ConstIterator( nbVertices() ); }
 
    
    //---------------------------------------------------------------------------
  public:
 
    const std::vector< RealPoint >& positions() const
    { return myPositions; }
 
    RealPoint& position( Vertex v )
    { return myPositions[ v ]; }
    
    const RealPoint& position( Vertex v ) const
    { return myPositions[ v ]; }
      
    const std::vector< RealVector >& vertexNormals() const
    { return myVertexNormals; }
 
    std::vector< RealVector >& vertexNormals() 
    { return myVertexNormals; }
 
    RealVector& vertexNormal( Vertex v )
    { return myVertexNormals[ v ]; }
    
    const RealVector& vertexNormal( Vertex v ) const
    { return myVertexNormals[ v ]; }
    
    const std::vector< RealVector >& faceNormals() const
    { return myFaceNormals; }
 
    std::vector< RealVector >& faceNormals() 
    { return myFaceNormals; }
 
    RealVector& faceNormal( Face f )
    { return myFaceNormals[ f ]; }
    
    const RealVector& faceNormal( Face f ) const
    { return myFaceNormals[ f ]; }
    
    Scalar averageEdgeLength() const;
 
    Scalar distance(const Vertex i, const Vertex j) const
    {
      //unrolling for compiler optimization
      const auto p=this->myPositions[ i ];
      const auto q=this->myPositions[ j ];
      return std::sqrt( (p[0]-q[0])*(p[0]-q[0]) + (p[1]-q[1])*(p[1]-q[1])+ (p[2]-q[2])*(p[2]-q[2]));
    }
    
    Scalar localWindow( Face f ) const;
 
    void perturbateWithUniformRandomNoise( Scalar p );
 
    void perturbateWithAdaptiveUniformRandomNoise( Scalar p );
 
    RealPoint edgeCentroid( Index e ) const;
 
    RealPoint faceCentroid( Index f ) const;
 
    Scalar faceArea( Index f ) const;
 
    Scalars getMaxWeights( Index v ) const;
    
    WeightedFaces
    computeFacesInclusionsInBall( Scalar r, Index f ) const;
    
    WeightedFaces
    computeFacesInclusionsInBall( Scalar r, Index f, RealPoint p ) const;
    
    std::tuple< Vertices, WeightedEdges, WeightedFaces >
    computeCellsInclusionsInBall( Scalar r, Index f ) const;
 
    std::tuple< Vertices, WeightedEdges, WeightedFaces >
    computeCellsInclusionsInBall( Scalar r, Index f, RealPoint p ) const;
    
    Scalar faceInclusionRatio( RealPoint p, Scalar r, Index f ) const;
 
    Scalar edgeInclusionRatio( RealPoint p, Scalar r, Index e ) const;
 
    Scalar vertexInclusionRatio( RealPoint p, Scalar r, Index v ) const;
 
 
    //---------------------------------------------------------------------------
  public:
 
    bool isFlippable( const Edge e ) const;
 
    VertexPair otherDiagonal( const Edge e ) const;
    
    void flip( const Edge e, bool recompute_face_normals = false );
    
 
    //---------------------------------------------------------------------------
  public:
    
    void computeNeighbors();
    void computeEdges();
 
    
    // ----------------------- Interface --------------------------------------
  public:
 
    void selfDisplay ( std::ostream & out ) const;
 
    bool isValid() const;
    
    // ------------------------- Protected Datas ------------------------------
  protected:
    std::vector< Vertices >     myIncidentVertices;
    std::vector< Faces >        myIncidentFaces;
    std::vector< RealPoint >    myPositions;
    std::vector< RealVector >   myVertexNormals;
    std::vector< RealVector >   myFaceNormals;
    std::vector< Faces >        myNeighborFaces;
    std::vector< Vertices >     myNeighborVertices;
    std::vector< VertexPair >   myEdgeVertices;
    std::map< VertexPair,Edge > myVertexPairEdge;
    std::vector< Faces >        myEdgeFaces;
    std::vector< Faces >        myEdgeRightFaces;
    std::vector< Faces >        myEdgeLeftFaces;
 
    // ------------------------- Private Datas --------------------------------
  private:
 
 
    // ------------------------- Internals ------------------------------------
  protected:
 
    void removeIndex( std::vector< Index >& v, Index i )
    {
      const std::size_t n = v.size();
      for ( std::size_t j = 0; j < n; j++ )
        if ( v[ j ] == i )
          {
            std::swap( v[ j ], v.back() );
            v.resize( n - 1 );
            return;
          }
      trace.error() << "[SurfaceMesh::removeIndex] Index " << i
                    << " is not in vector:";
      for ( auto e : v ) std::cerr << " " << e;
      std::cerr << std::endl;
    }
 
    void replaceIndex( std::vector< Index >& v, Index i, Index ri )
    {
      const std::size_t n = v.size();
      for ( std::size_t j = 0; j < n; j++ )
        if ( v[ j ] == i )
          {
            v[ j ] = ri;
            return;
          }
      trace.error() << "[SurfaceMesh::replaceIndex] Index " << i
                    << " (subs=" << ri << ") is not in vector:";
      for ( auto e : v ) std::cerr << " " << e;
      std::cerr << std::endl;
    }
 
    void addIndex( std::vector< Index >& v, Index i )
    {
      v.push_back( i );
    }
    
 
    static Scalar rand01()
    { return (Scalar) rand() / (Scalar) RAND_MAX; }
    
  }; // end of class SurfaceMesh
 
  template < typename TRealPoint, typename TRealVector >
  std::ostream&
  operator<< ( std::ostream & out,
               const SurfaceMesh<TRealPoint, TRealVector> & object );  
  
} // namespace DGtal
 
// Includes inline functions.
#include "SurfaceMesh.ih"
//                                                                           //
 
#endif // !defined SurfaceMesh_h
 
#undef SurfaceMesh_RECURSES
#endif // else defined(SurfaceMesh_RECURSES)