PConvexity.h

 
#pragma once
 
#if defined(PConvexity_RECURSES)
#error Recursive header files inclusion detected in PConvexity.h
#else // defined(PConvexity_RECURSES)
#define PConvexity_RECURSES
 
#if !defined PConvexity_h
#define PConvexity_h
 
// Inclusions
#include <iostream>
#include <vector>
#include "DGtal/base/Common.h"
#include "DGtal/kernel/CInteger.h"
#include "DGtal/kernel/CSpace.h"
#include "DGtal/geometry/volumes/ConvexityHelper.h"
 
namespace DGtal
{
  namespace detail
  {
    template < Dimension dim,
               typename TInteger = DGtal::int32_t >
    struct RecursivePConvexity {
      BOOST_CONCEPT_ASSERT(( concepts::CInteger<TInteger> ) );
      using Integer        = TInteger;
      using Point          = DGtal::PointVector< dim, Integer >;
      using ProjPoint      = DGtal::PointVector< dim-1, Integer >;
      using ProjPConvexity = DGtal::detail::RecursivePConvexity< dim - 1, Integer >;
 
      RecursivePConvexity( Dimension bd = dim )
      {
        init( bd );
      }
  
      void init( Dimension bd = dim )
      {
        for ( Dimension j = 0; j < bd; j++ )
          projp.push_back( ProjPConvexity( j ) );
      }
  
      static
      bool is0Convex( const std::vector< Point >& X, bool safe ) 
      {
        if ( X.empty() ) return true;
        // Build polytope according to internal integer type.
        if ( safe )
          {
            typedef typename DGtal::detail::ConvexityHelperInternalInteger< Integer, true >::Type
              InternalInteger;
            const auto P = ConvexityHelper< dim, Integer, InternalInteger >::
	      computeLatticePolytope( X, false, false );
            const std::size_t number_lattice_points_in_P = P.count();
            return number_lattice_points_in_P == X.size();
          }
        else
          {
            typedef typename DGtal::detail::ConvexityHelperInternalInteger< Integer, false >::Type
              InternalInteger;
            const auto P = ConvexityHelper< dim, Integer, InternalInteger >::
	      computeLatticePolytope( X, false, false );
            const std::size_t number_lattice_points_in_P = P.count();
            return number_lattice_points_in_P == X.size();
          }
      }
  
      bool isPConvex( const std::vector< Point >& X, bool safe ) const
      {
        if ( ! is0Convex( X, safe ) ) return false;
        for ( std::size_t j = 0; j < projp.size(); j++ )
          {
            const auto pi_j_X = project( X, j );
            if ( ! projp[ j ].isPConvex( pi_j_X, safe ) ) return false;
          }
        return true;
      }
 
      static
      double convexityMeasure( const std::vector< Point >& X, bool safe ) 
      {
        if ( X.empty() ) return 1.0;
        // Build polytope according to internal integer type.
        if ( safe )
          {
            typedef typename DGtal::detail::ConvexityHelperInternalInteger< Integer, true >::Type
              InternalInteger;
            const auto P = ConvexityHelper< dim, Integer, InternalInteger >::
	      computeLatticePolytope( X, false, false );
            const std::size_t number_lattice_points_in_P = P.count();
            return double( X.size() ) / double( number_lattice_points_in_P );
          }
        else
          {
            typedef typename DGtal::detail::ConvexityHelperInternalInteger< Integer, false >::Type
              InternalInteger;
            const auto P = ConvexityHelper< dim, Integer, InternalInteger >::
	      computeLatticePolytope( X, false, false );
            const std::size_t number_lattice_points_in_P = P.count();
            return double( X.size() ) / double( number_lattice_points_in_P );
          }
      }
      
      double fullConvexityMeasure( const std::vector< Point >& X, bool safe ) const
      {
        double m = convexityMeasure( X, safe );
        for ( std::size_t j = 0; j < projp.size(); j++ )
          {
            auto pX = project( X, j );
            m      *= projp[ j ].fullConvexityMeasure( pX, safe );
          }
        return m;
      }
  
      static
      ProjPoint project( const Point& p, Dimension a ) 
      {
        ProjPoint pp;
        Dimension j = 0;
        for ( Dimension i = 0; i < Point::dimension; i++ )
          if ( i != a ) pp[ j++ ] = p[ i ];
        return pp;
      }
 
      static
      std::vector< ProjPoint > project( const std::vector< Point >& p, Dimension a ) 
      {
        std::vector< ProjPoint > pp( p.size() );
        for ( std::size_t i = 0; i < p.size(); i++ )
          pp[ i ] = project( p[ i ], a );
        std::sort( pp.begin(), pp.end() );
        auto last = std::unique( pp.begin(), pp.end() );
        pp.erase( last, pp.end() );
        return pp;
      }
 
      std::vector< ProjPConvexity > projp;
    };
 
    template < typename TInteger >
    struct RecursivePConvexity< 1, TInteger> {
      BOOST_CONCEPT_ASSERT(( concepts::CInteger<TInteger> ) );
      using Integer     = TInteger;
      using Point       = PointVector< 1, Integer >;
 
      RecursivePConvexity( Dimension /* unused parameter in 1D specialization */ ) 
      {}
      
      static
      bool is0Convex( std::vector< Point > X, bool safe ) 
      {
        (void) safe;
        std::sort( X.begin(), X.end() );
        return X.empty()
          || ( ( Integer(X.back()[ 0 ]) - Integer(X.front()[ 0 ]) + Integer(1) )
               == Integer( X.size() ) );
      }
 
      static
      bool isPConvex( const std::vector< Point >& X, bool safe )
      {
        return is0Convex( X, safe );
      }
      
      static
      double convexityMeasure( std::vector< Point > X, bool safe )
      {
        (void) safe; //< not used in dimension 1.
        if ( X.empty() ) return 1.0;
        std::sort( X.begin(), X.end() );
        Integer nb =  Integer(X.back()[ 0 ]) - Integer(X.front()[ 0 ]) + Integer(1);
        return double( X.size() ) / double( nb );
      }
      
      static
      double fullConvexityMeasure( const std::vector< Point >& X, bool safe )
      {
        return convexityMeasure( X, safe );
      }
      
    };
 
  } // namespace detail 
  
    
  
  // template class PConvexity
  template < typename TSpace >
  class PConvexity
  {
    BOOST_CONCEPT_ASSERT(( concepts::CSpace< TSpace > ));
 
  public:
    typedef PConvexity<TSpace>       Self;
    typedef TSpace                   Space;
    typedef typename Space::Integer  Integer;
    typedef typename Space::Point    Point;
    
    static const Dimension dimension = Space::dimension;
    using RPConvexity = DGtal::detail::RecursivePConvexity< dimension, Integer >;
 
    // ----------------------- Standard services --------------------------------------
  public:
 
    ~PConvexity() = default;
 
    PConvexity( bool safe = false )
      : myRPC(), mySafe( safe )
    {}
 
 
    // ----------------------- Convexity services --------------------------------------
  public:
 
    bool is0Convex( const std::vector< Point >& X ) const
    {
      return myRPC.is0Convex( X, mySafe );
    }
 
    bool isPConvex( const std::vector< Point >& X ) const
    {
      return myRPC.isPConvex( X, mySafe );
    }
    
 
    // ----------------------- Measure services --------------------------------------
  public:
 
    double convexityMeasure( const std::vector< Point >& X ) const
    {
      return myRPC.convexityMeasure( X, mySafe );
    }
 
    double fullConvexityMeasure( const std::vector< Point >& X ) const
    {
      return myRPC.fullConvexityMeasure( X, mySafe );
    }
 
 
    // ----------------------- Interface --------------------------------------
  public:
 
    void selfDisplay ( std::ostream & out ) const
    {
      out << "[PConvexity dim=" << dimension
          << " safe=" << ( mySafe ? "True" : "False" )
          << " #bits<int>=" << ( sizeof( Integer ) * 8 ) << "]";
    }
 
    bool isValid() const
    {
      return dimension >= 1;
    }
 
 
    // ------------------------- Protected Datas ------------------------------
  protected:
 
    RPConvexity myRPC;
    
    bool mySafe;
    
    // ------------------------- Private Datas --------------------------------
  private:
 
    // ------------------------- Internals ------------------------------------
  private:
    
  }; // end of class PConvexity
 
 
  template <typename TKSpace>
  std::ostream&
  operator<< ( std::ostream & out,
               const PConvexity<TKSpace> & object )
  {
    object.selfDisplay( out );
    return out;
  }
 
 
} // namespace DGtal
 
 
// Includes inline functions.
 
//                                                                           //
 
#endif // !defined PConvexity_h
 
#undef PConvexity_RECURSES
#endif // else defined(PConvexity_RECURSES)