IntegralIntervals.h

 
#pragma once
#if defined(IntegralIntervals_RECURSES)
#error Recursive header files inclusion detected in IntegralIntervals.h
#else // defined(IntegralIntervals_RECURSES)
#define IntegralIntervals_RECURSES
 
#if !defined IntegralIntervals_h
#define IntegralIntervals_h
 
#include <vector>
#include <set>
#include "DGtal/base/Common.h"
#include "DGtal/kernel/CBoundedNumber.h"
 
namespace DGtal
{
 
  // template class IntegralIntervals
  template < typename TInteger >
  class IntegralIntervals
  {
  public:
    BOOST_CONCEPT_ASSERT(( concepts::CBoundedNumber< TInteger > ));
 
    typedef TInteger Integer;
    using Self      = IntegralIntervals< Integer >;
    using Interval  = std::pair<Integer,Integer>;
    using Container = std::vector< Interval >;
    using Size      = std::size_t;
    using CIterator = typename Container::iterator;
    using IntegerRange = std::vector< Integer >;
    
    IntegralIntervals() = default;
 
    IntegralIntervals( const Self & other ) = default;
 
    IntegralIntervals( Self&& other ) = default;
 
    Self& operator=( const Self & other ) = default;
 
    Self& operator=( Self&& other ) = default;
    
    template <typename InputIterator>
    IntegralIntervals( InputIterator it, InputIterator itE )
    {
      if ( it == itE ) return;
      Integer first = *it;
      Integer last  = *it;
      for ( ++it; it != itE; ++it )
        {
          Integer x = *it;
          if ( first <= x && x <= last ) continue;
          if ( x == last+1 )  { last = x;  continue; }
          if ( x == first-1 ) { first = x; continue; }
          insert( first, last );
          first = x;
          last  = x;
        }
      insert( first, last );
    }
 
    void clear() { myData.clear(); }
 
    Container& data() { return myData; }
    const Container& data() const { return myData; }
 
    bool empty() const { return myData.empty(); }
 
    Size size() const
    {
      Size nb = 0;
      for ( const auto& I : myData ) nb += 1 + I.second - I.first;
      return nb;
    }
 
    Size capacity() const
    {
      return myData.capacity();
    }
    
    bool isConvex() const
    {
      return myData.size() <= 1;
    }
    
    std::set<Integer> integerSet() const
    {
      std::set<Integer> S;
      for ( const auto& I : myData )
        for ( Integer x = I.first; x <= I.second; x++ )
          S.insert( x );
      return S;
    }
    std::vector<Integer> integerVector() const
    {
      std::vector<Integer> S;
      for ( const auto& I : myData )
        for ( Integer x = I.first; x <= I.second; x++ )
          S.push_back( x );
      return S;
    }
      
    Size count( Integer x ) const
    {
      if ( empty() ) return 0;
      Size i = 0;
      Size j = myData.size() - 1;
      while ( i <= j )
        {
          const Size m = (i+j)/2;
          const Interval& I = myData[ m ]; // I = [a,...,b]
          if ( x < I.first ) // x < a
            {
              if ( m == 0 ) return 0;
              j = m - 1;
            }
          else if ( I.second < x ) // b < x
            i = m + 1;
          else // a <= x <= b
            return 1;
        }
      return 0;
    }
 
    void insert( Integer i )
    {
      insert( Interval( i, i ) );
    }
    
    void insert( Integer f, Integer l )
    {
      insert( Interval( f, l ) );
    }
    
    void insert( const Interval& I )
    {
      // Search position of first element.
      auto it = lowerBound( I.first );
      if ( it == myData.end() ) // if you reach the end, just add the interval
        {
          myData.push_back( I );
          if ( myData.size() >= 2 ) extend( myData.end() - 2 );
        }
      else if ( I.first < it->first )
        {
          // See if interval must merge with previous
          if ( it != myData.begin() )
            {
              auto it_prev = it; --it_prev;
              if ( I.first <= it_prev->second + 1 )
                {
                  it_prev->second = I.second;
                  extend( it_prev );
                  return;
                }
            }
          Size idx = it - myData.begin();
          // std::cout << "(Inserting " << idx << ")" << std::endl;;
          myData.insert( it, I );
          extend( myData.begin() + idx );
        }
      else // it->first <= I.first <= it->second
        {
          it->second = std::max( it->second, I.second );
          extend( it );
        }
    }
 
    void erase( Integer i )
    {
      erase( Interval( i, i ) );
    }
    void erase( Integer f, Integer l )
    {
      erase( Interval( f, l ) );
    }
    
    void erase( const Interval& I )
    {
      for ( std::size_t i = 0; i < myData.size(); )
        {
          Interval& J = myData[ i ];
          // I=[a,b], J=[a',b'], a <= b, a' <= b'
          if ( I.second < J.first )
            { break; } // b < a' : no further intersection
          if ( J.second < I.first )
            { ++i; continue; } // b' < a : no further intersection
          // a' <= b and a <= b'
          //       a ----------  b
          // a' ...............  a'
          //       b' ................. b'
          //
          // a' ..................... b' => a'..a-1       b+1..b'
          Interval K1( J.first, I.first - 1 );
          Interval K2( I.second + 1, J.second );
          bool K1_exist = K1.second >= K1.first;
          bool K2_exist = K2.second >= K2.first;
          if ( K1_exist && K2_exist )
            {
              myData[ i ] = K2;
              myData.insert( myData.begin() + i, K1 );
              break; // no further intersection possible
            }
          else if ( K1_exist )
            {
              myData[ i ] = K1; i++;
            }
          else if ( K2_exist )
            {
              myData[ i ] = K2; break;
            }
          else
            {
              myData.erase( myData.begin() + i );
            }
        }
    }
 
    Self& add( const Self& other )
    {
      for ( const auto& I : other.myData )
        insert( I );
      return *this;
    }
    
    Self& subtract( const Self& other )
    {
      for ( const auto& I : other.myData )
        erase( I );
      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 );
    }
 
    Self starOfPoints() const
    {
      Self R( *this );
      for ( auto& I : R.myData )
        {
          I.first  = 2*I.first-1;
          I.second = 2*I.second+1;
        }
      return R;
    }
 
    Self starOfCells() const
    {
      Self R( *this );
      for ( auto i = 0; i < R.myData.size(); )
        {
          auto& I = R.myData[ i ];
          if ( ( I.first  & 0x1 ) == 0 ) I.first  -= 1;
          if ( ( I.second & 0x1 ) == 0 ) I.second += 1;
          // We have to be careful since extending this interval may
          // have reached the next interval.
          // We have to merge them in this case.
          i += 1;
          if ( i < R.myData.size() )
            {
              auto& Inext = R.myData[ i ];
              if ( Inext.first <= I.second+1 )
                {
                  I.second = Inext.second;
                  R.myData.erase( R.myData.begin() + i );
                  i -= 1;
                }
            }
        }
      return R;
    }
 
    IntegerRange extremaOfCells() const
    {
      IntegerRange C;
      for ( auto I : myData )
        {
          if ( ( I.first  & 0x1 ) != 0 ) I.first  -= 1;
          if ( ( I.second & 0x1 ) != 0 ) I.second += 1;
          for ( auto x = I.first; x <= I.second; x += 2 )
            C.push_back( x >> 1 ); // here x / 2 == x >> 1 since x is even
        }
      auto last = std::unique( C.begin(), C.end() );
      C.erase( last, C.end() );
      return C;
    }
    
    bool includes( const Self& other ) const
    {
      auto it = myData.cbegin();
      for ( const auto& I : other.myData )
        {
          // Find possible interval
          while ( it != myData.cend() && it->second < I.second ) ++it;
          if ( it == myData.cend() )  return false;
          if ( I.first < it->first ) return false;
        }
      return true;
    }
 
    bool equals( const Self& other ) const
    {
      if ( myData.size() != other.myData.size() ) return false;
      auto it = myData.cbegin();
      for ( const auto& I : other.myData )
        {
          if ( it->first != I.first || it->second != I.second ) return false;
          ++it;
        }
      return true;
    }
    
    // ----------------------- Interface --------------------------------------
public:
 
    void selfDisplay ( std::ostream & out ) const
    {
      out << "[";
      for ( const auto& I : myData )
        out << " (" << I.first << "," << I.second << ")";
      out << " ]";
    }
 
    bool isValid() const
    {
      for ( const auto& I : myData )
        if ( I.first > I.second ) return false;
      for ( Size i = 1; i < myData.size(); i++ )
        {
          if ( myData[i-1].second >= myData[i].first - 1 )
            return false;
        }
      return true;
    }
    
    // ------------------- protected services -------------------------------
  protected:
    
    void extend( CIterator it )
    {
      // std::cout  << "Extending" << std::endl;
      CIterator it_next = it; ++it_next;
      while ( it_next != myData.end() )
        {
          if ( it->second >= ( it_next->first - 1 ) )
            {
              it->second = std::max( it->second, it_next->second );
              ++it_next;
            }
          else break;
        }
      ++it;
      // std::cout  << "Erase from " << ( it - myData.begin() )
      //            << " to " << ( it_next - myData.begin() ) << std::endl;
      myData.erase( it, it_next );
    }
 
    CIterator lowerBound( Integer x ) 
    {
      // std::cout << "(lowerbound for " << x << ")" << std::endl;
      if ( empty() ) return myData.end();
      Size i = 0;
      Size j = myData.size() - 1;
      while ( i <= j )
        {
          const Size m = (i+j)/2;
          const Interval& I = myData[ m ]; // I = [a,...,b]
          if ( x < I.first ) // x < a
            {
              if ( m == 0 ) break;
              j = m - 1;
            }
          else if ( I.second < x ) // b < x
            i = m + 1;
          else // a <= x <= b
            return myData.begin() + m;
        }
      // std::cout << "(not found, return " << i <<  ")" << std::endl;
      return myData.begin() + i;
    }
    
    // ------------------- protected data -------------------------------
  protected:
    
    Container myData;
  };
 
  template <typename TInteger>
  std::ostream&
  operator<< ( std::ostream & out, const IntegralIntervals<TInteger> & object )
  {
    object.selfDisplay( out );
    return out;
  }
  
} // namespace DGtal
 
 
// Includes inline functions.
 
//                                                                           //
 
#endif // !defined IntegralIntervals_h
 
#undef IntegralIntervals_RECURSES
#endif // else defined(IntegralIntervals_RECURSES)