LatticeSetByIntervals.h

 
#pragma once
#if defined(LatticeSetByIntervals_RECURSES)
#error Recursive header files inclusion detected in LatticeSetByIntervals.h
#else // defined(LatticeSetByIntervals_RECURSES)
#define LatticeSetByIntervals_RECURSES
 
#if !defined LatticeSetByIntervals_h
#define LatticeSetByIntervals_h
 
#include <unordered_map>
#include <boost/iterator/iterator_facade.hpp>
#include <climits>
#include "DGtal/base/Common.h"
#include "DGtal/kernel/CUnsignedNumber.h"
#include "DGtal/kernel/CBoundedNumber.h"
#include "DGtal/kernel/CSpace.h"
#include "DGtal/kernel/PointVector.h"
#include "DGtal/kernel/IntegralIntervals.h"
 
namespace DGtal
{
 
  // template class LatticeSetByIntervals
  template < typename TSpace >
  class LatticeSetByIntervals
  {
  public:
    BOOST_CONCEPT_ASSERT(( concepts::CSpace< TSpace > ));
 
    typedef TSpace Space;
    using Self      = LatticeSetByIntervals< Space >;
    using Point     = typename Space::Point;
    using Vector    = typename Space::Vector;
    using Integer   = typename Space::Integer;
    using PointRange= std::vector< Point >;
    using Intervals = IntegralIntervals< Integer >;
    using Interval  = typename Intervals::Interval;
    using Container = std::map< Point, Intervals >;
    using RowIterator      = typename Container::iterator; 
    using RowConstIterator = typename Container::const_iterator;
    using LatticeSetByInterval = std::map< Point, Interval >;
    using Size      = std::size_t;
    using size_type = Size;
    static const Dimension dimension = Space::dimension;
 
    //------------------- standard services (construction, move) -------------------
  public:
    
    LatticeSetByIntervals( Dimension axis = 0 )
      : myAxis( axis ), myData() {}
 
    LatticeSetByIntervals( const Self & other ) = default;
 
    LatticeSetByIntervals( Self&& other ) = default;
 
    Self& operator=( const Self & other ) = default;
 
    Self& operator=( Self&& other ) = default;
    
    template <typename PointIterator>
    LatticeSetByIntervals( PointIterator it, PointIterator itE, Dimension axis = 0 )
      : myAxis( axis )
    {
      for ( ; it != itE; ++it ) {
        Point   q = *it;
        Integer x = q[ axis ];
        q[ axis ] = 0;
        myData[ q ].insert( x );
      }
    }
 
    LatticeSetByIntervals( const LatticeSetByInterval& aSet, Dimension axis )
      : myAxis( axis )
    {
      for ( const auto& aRow : aSet )
        myData[ aRow.first ].data().push_back( aRow.second );
    }
      
    void clear() { myData.clear(); }
    
    void setAxis( Dimension axis )
    {
      myData.clear();
      myAxis = axis;
    }
 
    Dimension axis() const
    { return myAxis; }
 
    Container& data() { return myData; }
 
 
    //------------------- conversion services -----------------------------
  public:
 
    PointRange toPointRange() const
    {
      PointRange X;
      X.reserve( size() );
      for ( const auto& pV : myData )
        {
          auto p    = pV.first;
          auto iVec = pV.second.integerVector();
          for ( auto x : iVec )
            {
              p[ myAxis ] = x;
              X.push_back( p );
            }
        }
      return X;
    }
    
    
    //------------------- capacity services -----------------------------
  public:
    
    bool empty() const
    {
      return myData.empty();
    }
      
    Size size() const
    {
      Size nb = 0;
      for ( const auto& pV : myData )
        nb += pV.second.size();
      return nb;
    }
      
    Size max_size() const
    {
      return myData.max_size();
    }
 
    
    //------------------- modifier services -----------------------------
  public:
 
    void insert( Point p )
    {
      Integer x   = p[ myAxis ];
      p[ myAxis ] = 0;
      myData[ p ].insert( x );
    }
 
    void erase( Point p )
    {
      Integer x   = p[ myAxis ];
      p[ myAxis ] = 0;
      auto it = myData.find( p );
      if ( it != myData.end() )
        {
          it->second.erase( x );
          if ( it->second.empty() )
            myData.erase( it ); // purge  element if it was the last in the row.
        }
    }
 
    
    void purge()
    {
      for ( auto it = myData.begin(), itE = myData.end(); it != itE; )
        if ( it->second.empty() )
          it = myData.erase( it );
        else ++it;
    }
    
 
    //------------------- set operations --------------------------------
  public:
 
    Self& add( const Self& other )
    {
      if ( other.axis() != axis() )
        {
          trace.error() << "[LatticeSetByInterval::add] "
                        << "Both lattice sets should share the same axis: "
                        << axis() << " != " << other.axis() << std::endl;
          return *this;
        }
      for ( const auto& pV : other.myData )
        {
          const Point& p = pV.first;
          auto it = myData.find( p );
          if ( it != myData.end() )
            it->second.add( pV.second );
          else
            myData[ p ] = pV.second;
        }
      return *this;
    }
    
    Self& subtract( const Self& other )
    {
      if ( other.axis() != axis() )
        {
          trace.error() << "[LatticeSetByInterval::subtract] "
                        << "Both lattice sets should share the same axis: "
                        << axis() << " != " << other.axis() << std::endl;
          return *this;
        }
      for ( const auto& pV : other.myData )
        {
          const Point& p = pV.first;
          auto it = myData.find( p );
          if ( it != myData.end() )
            it->second.subtract( pV.second );
        }
      return *this;
    }
 
    Self set_union( const Self& other ) const
    {
      Self U = *this;
      U.add( other );
      return U;
    }
 
    Self set_difference( const Self& other ) const
    {
      Self U = *this;
      U.subtract( other );
      return U;
    }
 
    Self set_intersection( const Self& other ) const
    {
      Self A_plus_B  = set_union( other );
      Self A_delta_B = set_symmetric_difference( other );
      return A_plus_B.subtract( A_delta_B );
    }
 
    Self set_symmetric_difference( const Self& other ) const
    {
      Self A_minus_B = *this;
      A_minus_B.subtract( other );
      Self B_minus_A = other;
      B_minus_A.subtract( *this );
      return A_minus_B.add( B_minus_A );
    }
 
    bool includes( const Self& other ) const
    {
      if ( other.axis() != axis() )
        {
          trace.error() << "[LatticeSetByInterval::subtract] "
                        << "Both lattice sets should share the same axis: "
                        << axis() << " != " << other.axis() << std::endl;
          return false;
        }
      for ( const auto& pV : other.myData )
        {
          const Point& p = pV.first;
          const auto  it = myData.find( p );
          if ( it == myData.cend() )                return false;
          if ( ! it->second.includes( pV.second ) ) return false;
        }
      return true;
    }
 
    bool equals( const Self& other ) const
    {
      if ( other.axis() != axis() )
        {
          trace.error() << "[LatticeSetByInterval::subtract] "
                        << "Both lattice sets should share the same axis: "
                        << axis() << " != " << other.axis() << std::endl;
          return false;
        }
      if ( myData.size() != other.myData.size() ) return false;
      auto it = myData.cbegin();
      for ( const auto& I : other.myData )
        {
          if ( it->first != I.first )            return false;
          if ( ! it->second.equals( I.second ) ) return false;
          ++it;
        }
      return true;
    }
 
    
    //------------------- topology operations --------------------------------
  public:
    
    Self starOfPoints() const
    {
      Self C( myAxis );
      // First step, place points as pointels and insert their star along
      // dimension a.
      for ( auto& pV : myData )
        {
          const Point q = 2 * pV.first;
          C.myData[ q ] = pV.second.starOfPoints();
        }
      // Second step, dilate along remaining directions
      for ( Dimension k = 0; k < dimension; k++ )
        {
          if ( k == myAxis ) continue;
          for ( const auto& value : C.myData )
            {
              Point    q = value.first;
              if ( q[ k ] & 0x1 ) continue;
              q[ k ]    -= 1;
              C.myData[ q ].add( value.second );
              q[ k ]    += 2;
              C.myData[ q ].add( value.second );
            }
        }
      return C;
    }
 
    Self starOfCells() const
    {
      Self C( *this );
      // First step, compute star along dimension a.
      for ( auto& pV : C.myData )
        pV.second = pV.second.starOfCells();
      // Second step, dilate along remaining directions
      for ( Dimension k = 0; k < dimension; k++ )
        {
          if ( k == myAxis ) continue;
          for ( const auto& value : C.myData )
            {
              Point    q = value.first;
              if ( q[ k ] & 0x1 ) continue;
              q[ k ]    -= 1;
              C.myData[ q ].add( value.second );
              q[ k ]    += 2;
              C.myData[ q ].add( value.second );
            }
        }
      return C;
    }
 
    Self skeletonOfCells() const
    {
      Self S( *this );
      // Now extracting implicitly its Skel
      for ( const auto& pV : myData )
        {
          const Point&   p = pV.first;
          const auto&    V = pV.second;
          for ( Dimension k = 0; k < dimension; k++ )
            {
              if ( k == myAxis ) continue;
              if ( ( p[ k ] & 0x1 ) != 0 ) continue; // if open along axis continue
              // if closed, check upper incident cells along direction k
              Point q = p;  q[ k ] -= 1;
              Point r = p;  r[ k ] += 1;
              auto itq = S.myData.find( q );
              if ( itq != S.myData.end() )
                {
                  auto& W = itq->second;
                  W.subtract( V );
                }
              auto itr = S.myData.find( r );
              if ( itr != S.myData.end() )
                {
                  auto& W = itr->second;
                  W.subtract( V );
                }
            }
        }
      // Extract skel along main axis
      for ( auto& value : S.myData )
        {
          auto & V = value.second;
          Intervals sub_to_V;
          for ( auto I : V.data() )
            {
              if ( ( I.first & 0x1 )  != 0 )
                {
                  if ( I.first != I.second )
                    sub_to_V.data().push_back( Interval{ I.first, I.first } );
                  I.first  += 1;
                }
              if ( ( I.second & 0x1 ) != 0 ) I.second -= 1;
              for ( auto x = I.first; x <= I.second; x += 2 )
                sub_to_V.data().push_back( Interval{ x+1, x+1 } );
            }
          V.subtract( sub_to_V );
        }
      // Erase empty stacks
      S.purge();
      return S;
    }
 
    PointRange extremaOfCells() const
    {
      typedef std::vector< Integer > Coordinates;
      std::map< Point, Coordinates > E;
      // Now extracting vertices along axis direction.
      for ( const auto& pV : myData )
        {
          const Point&   p = pV.first;
          const auto&    V = pV.second;
          E[ p ] = V.extremaOfCells();
        }
      std::map< Point, Coordinates > next_E;
      for ( Dimension k = 0; k != dimension; k++ )
        {
          if ( k == myAxis ) continue;
          for ( const auto& pC : E )
            {
              const auto& p = pC.first;
              Point       q = p;
              q[ k ]        = q[ k ] >> 1;
              bool odd      = ( p[ k ] & 0x1 ) != 0;
              { // odd/even always copy 
                auto it = next_E.find( q );
                if ( it == next_E.end() ) next_E[ q ] = pC.second;
                else
                  {
                    Coordinates F;
                    std::set_union( pC.second.cbegin(),  pC.second.cend(),
                                    it->second.cbegin(), it->second.cend(),
                                    std::back_inserter( F ) );
                    it->second = F;
                  }
              }
              if ( odd )
                { // odd: must copy forward also
                  q[ k ] += 1;
                  auto it = next_E.find( q );
                  if ( it == next_E.end() ) next_E[ q ] = pC.second;
                  else
                    {
                      Coordinates F;
                      std::set_union( pC.second.cbegin(),  pC.second.cend(),
                                      it->second.cbegin(), it->second.cend(),
                                      std::back_inserter( F ) );
                      it->second = F;
                    }
                }
            } // for ( const auto& pC : E )
          E.swap( next_E );
          next_E.clear();
        }
      // Build point range.
      PointRange R;
      for ( const auto& pC : E )
        {
          Point p = pC.first;
          for ( auto&& x : pC.second )
            {
              p[ myAxis ] = x;
              R.push_back( p );
            }
        }
      std::sort( R.begin(), R.end() );
      return R;
    }
    
    
    //------------------- specific services (interval insertion, removal) -------------
  public:
 
    size_type memory_usage() const noexcept
    {
      size_type nb = 0;
      for ( const auto& pV : myData )
        {
          nb += sizeof( Interval )  * pV.second.capacity() + sizeof( void* );
          nb += sizeof( pV ) + sizeof( void* );
        }
      nb += sizeof( Self );
      return nb;
    }
    
    Intervals& at( Point q )
    {
      q[ myAxis ] = 0;
      return myData[ q ];
    }
 
    void reset( Point p )
    {
      p[ myAxis ] = 0;
      auto it = myData.find( p );
      if ( it != myData.end() ) myData.erase( it );
    }
 
    
    Dimension myAxis;
    Container myData;
  };
 
} // namespace DGtal
 
 
// Includes inline functions.
 
//                                                                           //
 
#endif // !defined LatticeSetByIntervals_h
 
#undef LatticeSetByIntervals_RECURSES
#endif // else defined(LatticeSetByIntervals_RECURSES)