ConvexCellComplex.h

 
#pragma once
 
#if defined(ConvexCellComplex_RECURSES)
#error Recursive header files inclusion detected in ConvexCellComplex.h
#else // defined(ConvexCellComplex_RECURSES)
#define ConvexCellComplex_RECURSES
 
#if !defined ConvexCellComplex_h
#define ConvexCellComplex_h
 
// Inclusions
#include <iostream>
#include <string>
#include <vector>
#include "DGtal/base/Common.h"
#include "DGtal/kernel/PointVector.h"
#include "DGtal/math/linalg/SimpleMatrix.h"
#include "DGtal/geometry/volumes/BoundedLatticePolytope.h"
 
namespace DGtal
{
  // template class ConvexCellComplex
  template < typename TPoint >
  struct ConvexCellComplex {
    static const Dimension dimension = TPoint::dimension;
 
    typedef TPoint      Point;
    typedef std::size_t Index;
    typedef std::size_t Size;
 
    static const Index INFINITE_CELL = (Index) -1;
    typedef Index                    Cell;
    typedef std::pair< Index, bool > Face;
    typedef Index                    Vertex;
    typedef std::vector< Index >     IndexRange;
    typedef std::vector< Vertex >    VertexRange;
    typedef std::vector< Face >      FaceRange;
 
    typedef typename Point::Coordinate       Scalar;
    typedef PointVector< dimension, Scalar > Vector;
    typedef PointVector< dimension, double > RealPoint;
    typedef PointVector< dimension, double > RealVector;
    
    // ----------------------- standard services ---------------------------
  public:
 
    ConvexCellComplex()
    { clear(); }
 
    void clear()
    {
      cell_faces.clear();
      cell_vertices.clear();
      true_face_cell.clear();
      false_face_cell.clear();
      true_face_vertices.clear();
      vertex_position.clear();
      has_face_geometry = false;
      true_face_normal.clear();
      true_face_intercept.clear();
    }
    
    Size nbCells() const
    { return cell_faces.size(); }
    Size nbFaces() const
    { return true_face_cell.size(); }
    Size nbVertices() const
    { return vertex_position.size(); }
 
    // -------------------- primal structure topology services ---------------------
  public:
    
    Cell infiniteCell() const
    { return INFINITE_CELL; }
 
    bool isInfinite( const Cell c ) const
    { return c == INFINITE_CELL; }
      
    
    Face opposite( const Face f ) const
    { return std::make_pair( f.first, ! f.second ); }
 
    const FaceRange& cellFaces( const Cell c ) const
    { return cell_faces[ c ]; }
    
    const VertexRange& cellVertices( const Cell c ) const
    {
      ASSERT( c != INFINITE_CELL && c < nbCells() );
      ASSERT( ! cell_vertices.empty() );
      if ( cell_vertices.empty() )
        cell_vertices.resize( nbCells() );
      if ( cell_vertices[ c ].empty() )
        computeCellVertices( c );
      return cell_vertices[ c ];
    }
 
    const std::vector< VertexRange > & allCellVertices() const
    {
      if ( cell_vertices.empty() )
        cell_vertices.resize( nbCells() );
      for ( Index c = 0; c < nbCells(); ++c )
        if ( cell_vertices[ c ].empty() )
          computeCellVertices( c );
      return cell_vertices;
    }
    
    VertexRange faceVertices( const Face f ) const
    {
      if ( f.second ) return true_face_vertices[ f.first ];
      return VertexRange( true_face_vertices[ f.first ].crbegin(),
                          true_face_vertices[ f.first ].crend() );
    }
 
    Cell faceCell( const Face f ) const
    {
      return f.second ? true_face_cell[ f.first ] : false_face_cell[ f.first ];
    }
 
    VertexRange faceComplementVertices( const Face f ) const
    {
      VertexRange result;
      const Cell c = faceCell( f );
      if ( isInfinite( c ) ) return result;
      const auto& c_vtcs = cellVertices( c );
      auto        f_vtcs = true_face_vertices[ f.first ];
      std::sort( f_vtcs.begin(), f_vtcs.end() );
      for ( auto v : c_vtcs )
        if ( ! std::binary_search( f_vtcs.cbegin(), f_vtcs.cend(), v ) )
          result.push_back( v );
      return result;
    }
    
    // -------------------- primal structure geometry services ---------------------
  public:
 
    std::vector< Point > cellVertexPositions( const Cell c ) const
    {
      ASSERT( ! isInfinite( c ) );
      const auto vtcs = cellVertices( c );
      std::vector< Point > pts( vtcs.size() );
      std::transform( vtcs.cbegin(), vtcs.cend(), pts.begin(),
                      [&] ( Vertex v ) { return this->position( v ); } );
      return pts;
    }
 
    std::vector< Point > faceVertexPositions( const Face f ) const
    {
      const auto vtcs = faceVertices( f );
      std::vector< Point > pts( vtcs.size() );
      std::transform( vtcs.cbegin(), vtcs.cend(), pts.begin(),
                      [&] ( Vertex v ) { return this->position( v ); } );
      return pts;
    }
    
    Point position( const Vertex v ) const
    {
      return vertex_position[ v ];
    }
 
    RealPoint toReal( const Point p ) const
    {
      RealPoint x;
      for ( Dimension k = 0; k < dimension; ++k )
        x[ k ] = (double) p[ k ];
      return x;
    }
 
    template < typename LatticePoint >
    LatticePoint toLattice( const RealPoint p, double factor = 1.0 ) const
    {
      LatticePoint x;
      for ( Dimension k = 0; k < dimension; ++k )
        x[ k ] = (typename LatticePoint::Coordinate) round( p[ k ] * factor );
      return x;
    }
 
    
    RealPoint cellBarycenter( const Cell c ) const
    {
      ASSERT( ! isInfinite( c ) );
      RealPoint b;
      const auto& vtcs = cellVertices( c );
      for ( auto v : vtcs )
        b += toReal( position( v ) );
      return b / vtcs.size();
    }
 
    template < typename LatticePolytope >
    LatticePolytope
    cellLatticePolytope( const Cell c, const double factor = 1.0 ) const
    {
      ASSERT( hasFaceGeometry() );
      ASSERT( ! isInfinite( c ) );
      typedef typename LatticePolytope::Point     LatticePoint;
      typedef typename LatticePolytope::HalfSpace HalfSpace;
      typedef typename LatticePolytope::Domain    Domain;
      const auto vtcs = cellVertices( c );
      std::vector< LatticePoint > pts;
      for ( auto v : vtcs )
        pts.push_back( toLattice< LatticePoint >( position ( v ), factor ) );
      LatticePoint l = pts[ 0 ];
      LatticePoint u = pts[ 0 ];
      for ( const auto& p : pts ) {
        l = l.inf( p );
        u = u.sup( p );
      }
      Domain domain( l, u );
      std::vector< HalfSpace > HS;
      for ( auto f : cellFaces( c ) ) {
        HS.emplace_back( faceNormal( f ), faceIntercept( f ) );
      }
      return LatticePolytope( domain, HS.cbegin(), HS.cend(), false );
    }
 
    const Vector& faceNormal( const Face f ) const
    {
      ASSERT( hasFaceGeometry() );
      ASSERT( f.first < nbFaces() );
      if ( f.second ) return  true_face_normal[ f.first ];
      else            return -true_face_normal[ f.first ];
    }
 
    Scalar faceIntercept( const Face f ) const
    {
      ASSERT( hasFaceGeometry() );
      ASSERT( f.first < nbFaces() );
      if ( f.second ) return  true_face_intercept[ f.first ];
      else            return -true_face_intercept[ f.first ];
    }
      
    bool hasFaceGeometry() const
    { return has_face_geometry; } 
 
    void requireFaceGeometry() 
    {
      if ( ! hasFaceGeometry() )
        computeFaceGeometry();
    }
 
    void unrequireFaceGeometry()
    {
      has_face_geometry = false;
      true_face_normal.clear();
      true_face_intercept.clear();
    }
 
    
    // ----------------------- datas for primal structure ----------------------
  public:
    
    // for each cell, gives its faces
    std::vector< FaceRange >   cell_faces;
    // for each cell, gives its vertices
    mutable std::vector< VertexRange > cell_vertices;
    // for each 'true' face, gives its cell (or INFINITE_CELL)
    std::vector< Cell >        true_face_cell;
    // for each 'false' face, gives its cell (or INFINITE_CELL)
    std::vector< Cell >        false_face_cell;
    // for each true face, gives its vertices in order.
    std::vector< VertexRange > true_face_vertices;
    // for each vertex, gives its position
    std::vector< Point >       vertex_position;
 
    bool has_face_geometry;
    std::vector< Vector > true_face_normal;
    std::vector< Scalar > true_face_intercept;
    
    // ----------------------- Interface ----------------------
  public:
 
    void selfDisplay ( std::ostream & out ) const
    {
      out << "[ConvexCellComplex<" << dimension << ">"
          << " #C=" << nbCells()
          << " #F=" << nbFaces()
          << " #V=" << nbVertices();
      if ( hasFaceGeometry() ) out << " hasFaceGeometry";
      out << " ]";
    }
 
    bool isValid() const
    {
      return nbCells() != 0;
    }
 
    // ----------------------- protected services ----------------------
  protected:
 
    void computeCellVertices( const Cell c ) const
    {
      std::set< Vertex > vertices;
      for ( auto f : cell_faces[ c ] )
        vertices.insert( true_face_vertices[ f.first ].cbegin(),
                         true_face_vertices[ f.first ].cend() );
      cell_vertices[ c ] = VertexRange( vertices.cbegin(), vertices.cend() );
    }
 
    void computeFaceGeometry()
    {
      true_face_normal   .resize( nbFaces() );
      true_face_intercept.resize( nbFaces() );
      // Compute normals and intercepts
      for ( Index f = 0; f < nbFaces(); ++f ) {
        std::tie( true_face_normal[ f ], true_face_intercept[ f ] )
          = computeHalfSpace( true_face_vertices[ f ] );
        if ( true_face_normal[ f ] == Vector::zero )
          trace.error() << "[ConvexCellComplex::computeFaceGeometry]"
                        << " null normal vector at face " << f << std::endl;
      }
      // Orient them consistently
      for ( Index c = 0; c < nbCells(); c++ ) {
        //const auto& c_vtcs = cellVertices( c ); //not used
        for ( auto f : cellFaces( c ) ) {
          if ( ! f.second ) continue; // process only 'true' faces
          const   auto ov = faceComplementVertices( f );
          ASSERT( ! ov.empty() );
          const Vector  N = true_face_normal   [ f.first ];
          const Scalar nu = true_face_intercept[ f.first ];
          const Scalar iv = N.dot( position( ov.front() ) );
          if ( ( iv - nu ) > 0 ) {
            // Should be inside
            true_face_normal   [ f.first ] = -N;
            true_face_intercept[ f.first ] = -nu;
            std::reverse( true_face_vertices[ f.first ].begin(),
                          true_face_vertices[ f.first ].end() );
          }
          reorderFaceVertices( f.first );
        }
      }
      has_face_geometry = true;
    }
 
    std::pair< Vector, Scalar >
    computeHalfSpace( const VertexRange& v ) const
    {
      typedef DGtal::SimpleMatrix< Scalar, dimension, dimension > Matrix;
      Matrix A;
      Vector N;
      Scalar c;
      for ( int i = 1; i < dimension; i++ )
        for ( int j = 0; j < dimension; j++ )
          A.setComponent( i-1, j, position( v[ i ] )[ j ] - position( v[ 0 ] )[ j ] );
      for ( int j = 0; j < dimension; j++ )
        N[ j ] = A.cofactor( dimension-1, j );
      c = N.dot( position( v[ 0 ] ) );
      return std::make_pair( N, c );
    }
 
    void reorderFaceVertices( Index f )
    {
      VertexRange& result = true_face_vertices[ f ];
      Vector       N      = true_face_normal  [ f ];
      // Sort a facet such that its points, taken in order, have
      // always the same orientation of the facet.  More precisely,
      // facets span a `dimension-1` vector space. There are thus
      // dimension-2 fixed points, and the last ones (at least two)
      // may be reordered.
      VertexRange splx( dimension );
      for ( Dimension k = 0; k < dimension-2; k++ )
        splx[ k ] = result[ k ];
      std::sort( result.begin()+dimension-2, result.end(),
                 [&] ( Index i, Index j )
                 {
                   splx[ dimension-2 ] = i;
                   splx[ dimension-1 ] = j;
                   const auto H = computeHalfSpace( splx );
                   const auto orient = N.dot( H.first );
                   return orient > 0;
                 } );
    }
    
 
  };  // class ConvexCellComplex
 
  template <typename TPoint>
  std::ostream&
  operator<< ( std::ostream & out, const ConvexCellComplex<TPoint> & object )
  {
    object.selfDisplay( out );
    return out;
  }
  
} // namespace DGtal
 
#endif // !defined ConvexCellComplex_h
 
#undef ConvexCellComplex_RECURSES
#endif // else defined(ConvexCellComplex_RECURSES)