geometry/volumes/exampleBoundedLatticePolytopeCount4D.cpp

Computes the number of lattice points within random polytopes in two different ways: (1) by simple range scanning and inside test, (2) by computing intersections along an axis.

Second method is much faster (from 5x to easily 100x).

# 10 polytopes made from 10 points in range [-50:50]^4
./examples/geometry/volumes/exampleBoundedLatticePolytopeCount4D 10 10 50

outputs

New Block [Compute 4D polytopes]
EndBlock [Compute 4D polytopes] (1.68596 ms)
Computed 10 4D polytopes with 32.7 facets on average, in 1.68596 ms.
New Block [Compute number of lattice points within polytope (slow)]
EndBlock [Compute number of lattice points within polytope (slow)] (8301.91 ms)
New Block [Compute number of lattice points within polytope (fast)]
EndBlock [Compute number of lattice points within polytope (fast)] (597.177 ms)
Computed inside points is OK
Reference method computed 29885333 points in 8301.91 ms.
Fast method computed 29885333 points in 597.177 ms.

 
#include "DGtal/base/Common.h"
#include "DGtal/geometry/volumes/ConvexityHelper.h"
#include "DGtal/geometry/volumes/BoundedLatticePolytope.h"
 
int main( int argc, char* argv[] )
{
  int N  = argc > 1 ? atoi( argv[ 1 ] ) : 10;  // nb of polytopes
  int nb = argc > 2 ? atoi( argv[ 2 ] ) : 10;   // nb points per polytope
  int R  = argc > 3 ? atoi( argv[ 3 ] ) : 50;   // max diameter of shape
 
  typedef int64_t                              Integer;
  typedef DGtal::ConvexityHelper< 4, Integer > Helper;
  typedef Helper::LatticePolytope              Polytope;
  typedef Polytope::Point                      Point;
  // Compute all polytopes
  DGtal::trace.beginBlock( "Compute 4D polytopes" );
  std::vector< Polytope > polytopes;
  int sum_nb_facets = 0;
  for ( int i = 0; i < N; i++ )
    {
      std::vector< Point > V;
      for ( int i = 0; i < nb; i++ ) {
        Point p( rand() % (2*R+1) - R, rand() % (2*R+1) - R,
                 rand() % (2*R+1) - R, rand() % (2*R+1) - R );
        V.push_back( p );
      }
      Polytope P = Helper::computeLatticePolytope( V );
      sum_nb_facets += P.nbHalfSpaces();
      polytopes.push_back( P );
    }
  double t1 = DGtal::trace.endBlock();
  DGtal::trace.info() << "Computed " << N
               << " 4D polytopes with " << ( sum_nb_facets / (double) N )
               << " facets on average, in " << t1 << " ms." << std::endl;
  // Count interior points (slow method)
  DGtal::trace.beginBlock( "Compute number of lattice points within polytope (slow)" );
  std::size_t slow_nb = 0;
  std::vector< Integer > slow_counts;
  for ( const auto& P : polytopes )
    {
      const auto nb = P.countByScanning();
      slow_nb      += nb;
      slow_counts.push_back( nb ); 
    }
  double t2 = DGtal::trace.endBlock();
  // Count interior points (fast method)
  DGtal::trace.beginBlock( "Compute number of lattice points within polytope (fast)" );
  std::size_t fast_nb = 0;
  std::vector< Integer > fast_counts;
  for ( const auto& P : polytopes )
    {
      const auto nb = P.count();
      fast_nb      += nb;
      fast_counts.push_back( nb );
    }
  double t3 = DGtal::trace.endBlock();
  bool ok = std::equal( slow_counts.cbegin(), slow_counts.cend(), fast_counts.cbegin() );
  DGtal::trace.info() << "Computed inside points is " << ( ok ? "OK" : "ERROR" ) << std::endl;
  DGtal::trace.info() << "Reference method computed " << slow_nb
                      << " points in " << t2 << " ms." << std::endl;
  DGtal::trace.info() << "Fast method computed " << fast_nb
                      << " points in " << t3 << " ms." << std::endl;
  return 0;
}