NeighborhoodConvexityAnalyzer.h

 
#pragma once
 
#if defined(NeighborhoodConvexityAnalyzer_RECURSES)
#error Recursive header files inclusion detected in NeighborhoodConvexityAnalyzer.h
#else // defined(NeighborhoodConvexityAnalyzer_RECURSES)
#define NeighborhoodConvexityAnalyzer_RECURSES
 
#if !defined NeighborhoodConvexityAnalyzer_h
#define NeighborhoodConvexityAnalyzer_h
 
// Inclusions
#include <iostream>
#include <list>
#include <vector>
#include <string>
#include <bitset>
#include <unordered_set>
#include "DGtal/base/Common.h"
#include "DGtal/base/Clone.h"
#include "DGtal/base/ConstExpressions.h"
#include "DGtal/base/TimeStampMemoizer.h"
#include "DGtal/kernel/CPointPredicate.h"
#include "DGtal/kernel/CBoundedNumber.h"
#include "DGtal/kernel/domains/HyperRectDomain.h"
#include "DGtal/topology/CCellularGridSpaceND.h"
#include "DGtal/topology/KhalimskySpaceND.h"
#include "DGtal/geometry/volumes/DigitalConvexity.h"
 
namespace DGtal
{
  
  // template class NeighborhoodConvexityAnalyzer
  template < typename TKSpace, int K >
  class NeighborhoodConvexityAnalyzer
  {
    BOOST_CONCEPT_ASSERT(( concepts::CCellularGridSpaceND< TKSpace > ));
 
  public:
    typedef NeighborhoodConvexityAnalyzer<TKSpace,K> Self;
    typedef TKSpace                         KSpace;
    typedef typename KSpace::Space          Space;
    typedef typename KSpace::Integer        Integer;
    typedef typename KSpace::Point          Point;
    typedef typename KSpace::Vector         Vector;
    typedef typename KSpace::Cell           Cell;
    typedef std::vector<Point>              PointRange;
    typedef HyperRectDomain<Space>          Domain;
    typedef std::size_t                     Size;
    
    static const Dimension dimension  = KSpace::dimension;
    static const Size      neigh_size = functions::const_pow( 2*K+1, dimension ); 
    static const Size      middle     = functions::const_middle( K, dimension );
    static const bool  false_positive = ( dimension > 2 ) || ( K > 1 );
 
    typedef std::bitset< functions::const_pow( 2*K+1, dimension ) > Configuration;
    typedef std::bitset< 9 > BasicConfiguration;
 
    
    // ------------------------- Standard services --------------------------------
  public:
 
    ~NeighborhoodConvexityAnalyzer() = default;
 
    NeighborhoodConvexityAnalyzer() = default;
    
    NeighborhoodConvexityAnalyzer ( const Self & other ) = default;
 
    NeighborhoodConvexityAnalyzer( Clone<KSpace> aKSpace, Size memoizer_size = 0 )
      : myDigConv( aKSpace ), myMemoizer( memoizer_size )
    {
      myDomain       = Domain( myDigConv.space().lowerBound(),
                               myDigConv.space().upperBound() );
      myComputations = 0;
      myResults      = 0;
      computeBasicFullConvexityTable();
      trace.info() << "Size=" << size() << " middle=" << middle << std::endl;
    }
 
    NeighborhoodConvexityAnalyzer( Point lo, Point hi,
                                   Size memoizer_size = 0 )
      : myDomain( lo, hi ), myDigConv( lo, hi ), myMemoizer( memoizer_size )
    {
      myComputations = 0;
      myResults      = 0;
      computeBasicFullConvexityTable();
      trace.info() << "Size=" << size() << " middle=" << middle << std::endl;
    }
 
    Self & operator= ( const Self & other ) = default;
 
    const KSpace& space() const
    {
      return myDigConv.space();
    }
 
    const Domain& domain() const
    {
      return myDomain;
    }
 
    static int radius()
    { return K; }
 
    static Size size()
    { return neigh_size; }
    
 
    // -------------------- Neighborhood and convexity  services -----------------------
  public:
 
    template < typename PointPredicate >
    void setCenter( Point c, const PointPredicate& X );
 
    Point center() const
    {
      return myCenter;
    }
    
    Configuration configuration() const
    {
      return myCfgX;
    }
 
    bool isCenterInX() const
    {
      return myCenterInX;
    }
 
    bool isFullyConvexCollapsible()
    {
      if ( isCenterInX() )
        return ( myNbInX >= 1 ) // ( ! myLocalX.empty() )
          && isFullyConvex( true )
          && isFullyConvex( false );
      else
        return ( size() - myNbInX >= 2 ) // ( ! myLocalCompX.empty() )
          && isComplementaryFullyConvex( true )
          && isComplementaryFullyConvex( false );
    }
 
    bool isLikelyNoise()
    {
      if ( isCenterInX() )
        return ( myNbInX >= 1 ) // ( ! myLocalX.empty() )
          && ! isFullyConvex( true )
          && isFullyConvex( false )
          && isComplementaryFullyConvex( true );
      else
        return ( size() - myNbInX >= 2 ) // ( ! myLocalCompX.empty() )
          && ! isComplementaryFullyConvex( true )
          && isComplementaryFullyConvex( false )
          && isFullyConvex( true );
    }
    
    bool is0ConvexCollapsible()
    {
      if ( isCenterInX() )
        return ( myNbInX >= 1 ) // ( ! myLocalX.empty() )
          && is0Convex( true )
          && is0Convex( false );
      else
        return ( size() - myNbInX >= 2 ) // ( ! myLocalCompX.empty() )
          && isComplementary0Convex( true )
          && isComplementary0Convex( false );
    }
 
    bool isFullyConvex( bool with_center )
    {
      int mask = with_center
        ? FullConvexity_X_with_center : FullConvexity_X_without_center;
      if ( myComputations & mask ) return bool( myResults & mask );
      bool ok;
      // Check memoizer
      if ( myMemoizer.isValid() )
        {
          auto cfg = makeConfiguration( myCfgX, false, with_center );
          auto   p = myMemoizer.get( cfg );
          ok       = p.first; // may not be correct
          bool memoized = p.second;
          if ( ! memoized )
            {
              // Need to compute full convexity property
              ok = checkBasicConfigurationsFullConvexity( false, with_center );
              if ( ok && false_positive )
                { // need to do the true computation.
                  std::vector< Point > localX;
                  getLocalX( localX, with_center );
                  ok = myDigConv.isFullyConvex( localX );
                }
              myMemoizer.set( cfg, ok );
            }
        }
      else
        {
          ok = checkBasicConfigurationsFullConvexity( false, with_center );
          if ( ok && false_positive )
            { // need to do the true computation.
              std::vector< Point > localX;
              getLocalX( localX, with_center );
              ok = myDigConv.isFullyConvex( localX );
            }
        }
      // auto cfg = makeConfiguration( myCfgX, false, with_center );
      // // Check memoizer
      // if ( myMemoizer.isValid() )
      //   {
      //     auto   p = myMemoizer.get( cfg );
      //     ok       = p.first; // may not be correct
      //     memoized = p.second;
      //   }
      // if ( ! memoized )
      //   {
      //     // Need to compute full convexity property
      //     ok = checkBasicConfigurationsFullConvexity( false, with_center );
      //     if ( ok && false_positive )
      //       { // need to do the true computation.
      //         std::vector< Point > localX;
      //         getLocalX( localX, with_center );
      //         ok = myDigConv.isFullyConvex( localX );
      //       }
      //     if ( myMemoizer.isValid() )
      //       myMemoizer.set( cfg, ok );
      //   }
      myComputations |= mask;
      if ( ok ) myResults |= mask;
      return ok;
    }
 
    bool isComplementaryFullyConvex( bool with_center )
    {
      int mask = with_center
        ? FullConvexity_CompX_with_center : FullConvexity_CompX_without_center;
      if ( myComputations & mask ) return bool( myResults & mask );
      bool ok;
      // Check memoizer
      if ( myMemoizer.isValid() )
        {
          auto cfg = makeConfiguration( myCfgX, true, with_center );
          auto   p = myMemoizer.get( cfg );
          ok       = p.first; // may not be correct
          bool memoized = p.second;
          if ( ! memoized )
            {
              // Need to compute full convexity property
              ok = checkBasicConfigurationsFullConvexity( true, with_center );
              if ( ok && false_positive )
                { // need to do the true computation.
                  std::vector< Point > localCompX;
                  getLocalCompX( localCompX, with_center );
                  ok = myDigConv.isFullyConvex( localCompX );
                }
              myMemoizer.set( cfg, ok );
            }
        }
      else
        {
          // Need to compute full convexity property
          ok = checkBasicConfigurationsFullConvexity( true, with_center );
          if ( ok && false_positive )
            { // need to do the true computation.
              std::vector< Point > localCompX;
              getLocalCompX( localCompX, with_center );
              ok = myDigConv.isFullyConvex( localCompX );
            }
        }
      // bool memoized = false;
      // auto cfg = makeConfiguration( myCfgX, true, with_center );
      // // Check memoizer
      // if ( myMemoizer.isValid() )
      //   {
      //     auto   p = myMemoizer.get( cfg );
      //     ok       = p.first; // may not be correct
      //     memoized = p.second;
      //   }
      // if ( ! memoized )
      //   {
      //     // Need to compute full convexity property
      //     ok = checkBasicConfigurationsFullConvexity( true, with_center );
      //     if ( ok && false_positive )
      //       { // need to do the true computation.
      //         std::vector< Point > localCompX;
      //         getLocalCompX( localCompX, with_center );
      //         ok = myDigConv.isFullyConvex( localCompX );
      //       }
      //     if ( myMemoizer.isValid() )
      //       myMemoizer.set( cfg, ok );
      //   }
      myComputations |= mask;
      if ( ok ) myResults |= mask;
      return ok;
    }
 
    bool is0Convex( bool with_center )
    {
      int mask = with_center
        ? FullConvexity_X_with_center : FullConvexity_X_without_center;
      if ( myComputations & mask ) return bool( myResults & mask );
      // Need to compute full convexity property
      bool ok = checkBasicConfigurations0Convexity( false, with_center );
      if ( ok && false_positive )
        { // need to do the true computation.
          std::vector< Point > localX;
          getLocalX( localX, with_center );
          ok = myDigConv.is0Convex( localX );
        }
      myComputations |= mask;
      if ( ok ) myResults |= mask;
      return ok;
    }
 
    bool isComplementary0Convex( bool with_center )
    {
      int mask = with_center
        ? FullConvexity_CompX_with_center : FullConvexity_CompX_without_center;
      if ( myComputations & mask ) return bool( myResults & mask );
      // Need to compute full convexity property
      bool ok = checkBasicConfigurations0Convexity( true, with_center );
      if ( ok && false_positive )
        { // need to do the true computation.
          std::vector< Point > localCompX;
          getLocalCompX( localCompX, with_center );
          ok = myDigConv.is0Convex( localCompX );
        }
      myComputations |= mask;
      if ( ok ) myResults |= mask;
      return ok;
    }
 
    static Configuration makeConfiguration( Configuration current,
                                            bool complement, bool with_center )
    {
      if ( complement )
        {
          current.flip(); // current = ~current;
          if ( ! with_center ) current.reset( middle );
          else current.set( middle );
        }
      else
        {
          if ( with_center ) current.set( middle );
        }
      return current;
    }
    
    void getLocalX( std::vector< Point >& localX, bool with_center ) const;
 
    void getLocalCompX( std::vector< Point >& localCompX, bool with_center ) const;
    
    
    // ------------------------- Protected Datas ------------------------------
  protected:
    Domain myDomain;
    DigitalConvexity<KSpace> myDigConv;
    Point myCenter;
    TimeStampMemoizer< Configuration, bool > myMemoizer;
    bool myCenterInX;
    Size myNbInX;
    Configuration myCfgX;
    std::vector< BasicConfiguration > myBasicCfgX;
 
    std::bitset< 512 > myBasicFullConvexityTable;
    std::bitset< 512 > myBasic0ConvexityTable;
    
    enum Computation {
      FullConvexity_X_with_center        = 0x1,
      FullConvexity_X_without_center     = 0x2,
      FullConvexity_CompX_with_center    = 0x4,
      FullConvexity_CompX_without_center = 0x8,
      Convexity_X_with_center            = 0x10,
      Convexity_X_without_center         = 0x20,
      Convexity_CompX_with_center        = 0x40,
      Convexity_CompX_without_center     = 0x80,
    };
    int myComputations;
    int myResults;
      
    // ------------------------- Private Datas --------------------------------
  private:
 
    // ------------------------- Internals ------------------------------------
  private:
    void computeBasicFullConvexityTable();
 
    bool checkBasicConfigurationsFullConvexity
    ( bool compX, bool with_center ) const;
 
    bool checkBasicConfigurations0Convexity
    ( bool compX, bool with_center ) const;
 
    void computeBasicConfigurations
    ( Configuration cfg, std::vector< BasicConfiguration > & result ) const;
 
    BasicConfiguration computeCentralBasicConfiguration
    ( Configuration cfg, Dimension i, Dimension j ) const;
 
    
  }; // end of class NeighborhoodConvexityAnalyzer
 
 
} // namespace DGtal
 
 
// Includes inline functions.
#include "NeighborhoodConvexityAnalyzer.ih"
 
//                                                                           //
 
#endif // !defined NeighborhoodConvexityAnalyzer_h
 
#undef NeighborhoodConvexityAnalyzer_RECURSES
#endif // else defined(NeighborhoodConvexityAnalyzer_RECURSES)