DGtal  1.4.beta
geometry/volumes/exampleBoundedLatticePolytopeCount3D.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).

# 100 polytopes made from 10 points in range [-50:50]^3
./examples/geometry/tools/exampleBoundedLatticePolytopeCount3D 100 10 50

outputs

New Block [Compute polytope]
EndBlock [Compute polytope] (14.6467 ms)
Computed 100 polytopes with 19.11 facets on average, in 14.6467 ms.
New Block [Compute number of lattice points within polytope (slow)]
EndBlock [Compute number of lattice points within polytope (slow)] (853.593 ms)
New Block [Compute number of lattice points within polytope (fast)]
EndBlock [Compute number of lattice points within polytope (fast)] (61.2751 ms)
Computed inside points is OK
Reference method computed 15234308 points in 853.593 ms.
Fast method computed 15234308 points in 61.2751 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 ] ) : 100; // 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 Helper::LatticePolytope Polytope;
// Compute all polytopes
DGtal::trace.beginBlock( "Compute 3D 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 );
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
<< " 3D 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;
}
void beginBlock(const std::string &keyword="")
std::ostream & info()
double endBlock()
boost::int64_t int64_t
signed 94-bit integer.
Definition: BasicTypes.h:74
Trace trace
Definition: Common.h:153
Aim: Provides a set of functions to facilitate the computation of convex hulls and polytopes,...
int main(int argc, char **argv)
MyPointD Point
Definition: testClone2.cpp:383