SymmetricConvexExpander.h

 
#pragma once
 
#if defined(SymmetricConvexExpander_RECURSES)
#error Recursive header files inclusion detected in SymmetricConvexExpander.h
#else // defined(SymmetricConvexExpander_RECURSES)
#define SymmetricConvexExpander_RECURSES
 
#if !defined SymmetricConvexExpander_h
#define SymmetricConvexExpander_h
 
// Inclusions
#include <iostream>
#include "DGtal/base/Common.h"
#include "DGtal/topology/CCellularGridSpaceND.h"
#include "DGtal/geometry/volumes/DigitalConvexity.h"
 
namespace DGtal
{
 
  // class SymmetricConvexExpander
  template < typename TKSpace,
             typename TPointPredicate >
  class SymmetricConvexExpander
  {
    BOOST_CONCEPT_ASSERT(( concepts::CCellularGridSpaceND< TKSpace > ));
 
  public:
    typedef DigitalConvexity<TKSpace>       Self;
    typedef TKSpace                         KSpace;
    typedef TPointPredicate                 PointPredicate;
    typedef typename KSpace::Integer        Integer;
    typedef typename KSpace::Point          Point;
    typedef typename KSpace::Vector         Vector;
    typedef typename KSpace::Space          Space;
    typedef std::size_t                     Size;
    typedef DGtal::BoundedLatticePolytope < Space > LatticePolytope;
    typedef DGtal::BoundedRationalPolytope< Space > RationalPolytope;
    typedef std::vector<Point>              PointRange;
    typedef std::vector<Vector>             VectorRange;
    typedef std::unordered_set<Point>       PointSet;
    
    static const Dimension dimension = KSpace::dimension;
 
    typedef std::pair< Point, Integer >     Node;
 
    // Inversion order since priority queue output max element.
    struct NodeComparator {
      NodeComparator() = default;
      // p < q iff p.second < q.second
      bool operator()( const Node& p, const Node& q ) const
      {
        return p.second > q.second;
      }
    };
 
    typedef std::priority_queue< Node, std::vector<Node>, NodeComparator > NodeQueue;
 
    SymmetricConvexExpander( const PointPredicate& predicate,
                             const Point& kcenter,
                             const Point& lo, 
                             const Point& hi )
      : myPredicate( &predicate ), myConvexity( lo, hi )
    {
      init( kcenter );
    }
 
    bool predicate( const Point& p ) const
    {
      return (*myPredicate)( p );
    }
    
    void init( const Point& kcenter )
    {
      myKCenter = kcenter;
      myPoints.clear();
      myQ = NodeQueue();
      myM.clear();
      myPerfectSymmetry       = true;
      myPerfectSymmetryRadius = 0;
      // The starting points depend on the parity of the coordinates of the center.
      // There are from 1 to 2^d starting points.
      PointRange points;
      const auto x = myKCenter[ 0 ];
      if ( x % 2 == 0 )
        points.push_back( Point::base( 0, x / 2 ) );
      else
        {
          points.push_back( Point::base( 0, (x-1) / 2 ) );
          points.push_back( Point::base( 0, (x+1) / 2 ) );
        }
      for ( Dimension k = 1; k < dimension; k++ )
        {
          const auto n = points.size();
          const auto y = myKCenter[ k ];
          if ( y % 2 == 0 )
            {
              for ( auto i = 0; i < n; i++ )
                points[ i ][ k ] = y / 2;
            }
          else
            {
              points.resize( 2*n );
              const auto z  = (y-1)/2;
              const auto z1 = z + 1;
              for ( auto i = 0; i < n; i++ )
                {
                  points[ i ][ k ] = z;
                  Point q = points[ i ];
                  q[ k ]  = z1;
                  points[ i+n ]    = q;
                }                  
            }
        }
      // Keep only the points that satisfy the predicate.
      for ( auto&& p : points )
        {
          const Point sp = symmetric( p );
          if ( ! myM.count( p )
               && predicate( p ) && predicate( sp ) )
            {
              Node n( p, (2*p - myKCenter).squaredNorm() );
              myQ.push( n );
              myM.insert( p );
              myM.insert( sp );
            }
        }
    }
 
    bool advance( bool enforce_full_convexity )
    {
      while ( ! finished() )
        {
          const auto p = current().first;
          //NOT USED const auto d = current().second; // current ring distance
          const auto sp = symmetric( p );
          myPoints.insert( p );
          myPoints.insert( sp );
          PointRange X( myPoints.cbegin(), myPoints.cend() );
          if ( enforce_full_convexity && ! myConvexity.isFullyConvex( X ) )
            {
              myPoints.erase( p );
              myPoints.erase( sp );
              ignore();
            }
          else
            {
              expand();
              return true;
            }
        }
      return false;
    }
    
    const Node& current() const
    {
      return myQ.top();
    }
 
    void ignore()
    {
      myQ.pop();
    }
 
    void expand()
    {
      const Point p = current().first;
      myQ.pop();
      myPoints.insert( p );
      myPoints.insert( symmetric( p ) );
      const auto next_points = next( p );
      for ( auto&& q : next_points )
        {
          if ( ! myM.count( q ) )
            {
              const auto  sq   = symmetric( q );
              const bool  q_in = predicate( q );
              const bool sq_in = predicate( sq );
              if ( q_in && sq_in )
                {
                  Node n( q, (2*q - myKCenter).squaredNorm() );
                  myQ.push( n );
                  myM.insert(  q );
                  myM.insert( sq );
                  if ( myPerfectSymmetry )
                    myPerfectSymmetryRadius = std::max( myPerfectSymmetryRadius,
                                                        n.second );
                }
              else if ( ( q_in && ! sq_in ) || ( ! q_in && sq_in ) )
                {
                  myPerfectSymmetry = false;
                }
            }
        }
    }
 
    bool finished() const
    {
      return myQ.empty();
    }
    
    
    Point symmetric( const Point& p ) const
    {
      return myKCenter - p;
    }
 
    PointRange next( const Point& p ) const
    {
      PointRange N;
      Point d = 2*p - myKCenter;
      for ( Dimension i = 0; i < dimension; i++ )
        {
          if ( d[ i ] >= 0 ) N.push_back( p + Point::base( i, 1 ) );
          if ( d[ i ] <= 0 ) N.push_back( p - Point::base( i, 1 ) );
        }
      return N;
    }
 
    const PointPredicate* myPredicate;
 
    DigitalConvexity< KSpace > myConvexity;
    
    Point myKCenter;
 
    PointSet myPoints;
 
    NodeQueue  myQ;
    PointSet   myM;
 
    bool myPerfectSymmetry;
    Integer myPerfectSymmetryRadius;
    
  };
 
  
} // namespace DGtal
 
 
//                                                                           //
 
#endif // !defined SymmetricConvexExpander_h
 
#endif // !defined SymmetricConvexExpander_RECURSES