This example compares the speed of computation of P-convexity wrt to the computation of full convexity. Both definitions are equivalent but P-convexity is faster to compute, especially in higher dimensions.
Simply run the benchmark (it will take more than 1 hour on a M2 pro chip). It produces 9 files "timings-p-convexity-Z[d].txt", "timings-fc-convexity-Z[d].txt", and "timings-fcf-convexity-Z[d].txt", corresponding to P-convexity/full convexity/fast full convexity computation in Z[d]. Each data is a triplet (number of points, timings in ms, isConvex).
#include <iostream>
#include <vector>
#include <algorithm>
#include <chrono>
#include "DGtal/base/Common.h"
#include "DGtal/kernel/SpaceND.h"
#include "DGtal/kernel/domains/HyperRectDomain.h"
#include "DGtal/kernel/sets/DigitalSetBySTLSet.h"
#include "DGtal/topology/KhalimskySpaceND.h"
#include "DGtal/shapes/Shapes.h"
#include "DGtal/geometry/volumes/PConvexity.h"
#include "DGtal/geometry/volumes/DigitalConvexity.h"
using namespace std;
double rand01() {
return double( rand() ) / double( RAND_MAX ); }
template <Dimension dim>
void
std::size_t nb_tries, std::size_t nb_vertices, std::size_t range,
double pconvexity_probability = 0.5 )
{
typedef KhalimskySpaceND<dim,int64_t>
KSpace;
typedef PConvexity< Space > PConvexity;
DConvexity dconv( Point::diagonal( -1 ), Point::diagonal( range ) );
PConvexity pconv;
Domain domain( Point::diagonal( 0 ), Point::diagonal( range ) );
std::cout <<
"Computing " << nb_tries <<
" P-convexities in Z" <<
dim << std::endl;
for ( auto n = 0; n < nb_tries; ++n )
{
std::vector< Point > V;
for ( auto i = 0; i < nb_vertices; i++ ) {
for (
auto j = 0; j <
dim; j++ ) p[ j ] = rand() % range;
V.push_back( p );
}
std::vector< Point > X;
bool force_pconvexity =
rand01() < pconvexity_probability;
if ( force_pconvexity )
else
{
auto P = dconv.
CvxH( V );
}
std::chrono::high_resolution_clock::time_point
t1 = std::chrono::high_resolution_clock::now();
bool is_pconvex = pconv.isPConvex( X );
std::chrono::high_resolution_clock::time_point
t2 = std::chrono::high_resolution_clock::now();
double dt = std::chrono::duration_cast<std::chrono::nanoseconds>(t2 - t1).count();
results.push_back( std::make_tuple( X.size(),
dt/1e6, is_pconvex ) );
if ( force_pconvexity && ! is_pconvex )
trace.
warning() <<
"Invalid computation of either FC* or P-convexity !" << std::endl;
}
}
template <Dimension dim>
void
std::size_t nb_tries, std::size_t nb_vertices, std::size_t range,
double fconvexity_probability = 0.5 )
{
typedef KhalimskySpaceND<dim,int64_t>
KSpace;
typedef PConvexity< Space > PConvexity;
DConvexity dconv( Point::diagonal( -1 ), Point::diagonal( range ) );
PConvexity pconv;
Domain domain( Point::diagonal( 0 ), Point::diagonal( range ) );
std::cout <<
"Computing " << nb_tries <<
" full convexities in Z" <<
dim << std::endl;
for ( auto n = 0; n < nb_tries; ++n )
{
std::vector< Point > V;
for ( auto i = 0; i < nb_vertices; i++ ) {
for (
auto j = 0; j <
dim; j++ ) p[ j ] = rand() % range;
V.push_back( p );
}
std::vector< Point > X;
bool force_fconvexity =
rand01() < fconvexity_probability;
if ( force_fconvexity )
else
{
auto P = dconv.
CvxH( V );
}
std::chrono::high_resolution_clock::time_point
t1 = std::chrono::high_resolution_clock::now();
std::chrono::high_resolution_clock::time_point
t2 = std::chrono::high_resolution_clock::now();
double dt = std::chrono::duration_cast<std::chrono::nanoseconds>(t2 - t1).count();
results.push_back( std::make_tuple( X.size(),
dt/1e6, is_fconvex ) );
if ( force_fconvexity && ! is_fconvex )
trace.
warning() <<
"Invalid computation of either FC* or full convexity !" << std::endl;
}
}
template <Dimension dim>
void
std::size_t nb_tries, std::size_t nb_vertices, std::size_t range,
double fconvexity_probability = 0.5 )
{
typedef KhalimskySpaceND<dim,int64_t>
KSpace;
typedef PConvexity< Space > PConvexity;
DConvexity dconv( Point::diagonal( -1 ), Point::diagonal( range ) );
PConvexity pconv;
Domain domain( Point::diagonal( 0 ), Point::diagonal( range ) );
std::cout <<
"Computing " << nb_tries <<
" full convexities in Z" <<
dim << std::endl;
for ( auto n = 0; n < nb_tries; ++n )
{
std::vector< Point > V;
for ( auto i = 0; i < nb_vertices; i++ ) {
for (
auto j = 0; j <
dim; j++ ) p[ j ] = rand() % range;
V.push_back( p );
}
std::vector< Point > X;
bool force_fconvexity =
rand01() < fconvexity_probability;
if ( force_fconvexity )
else
{
auto P = dconv.
CvxH( V );
}
std::chrono::high_resolution_clock::time_point
t1 = std::chrono::high_resolution_clock::now();
std::chrono::high_resolution_clock::time_point
t2 = std::chrono::high_resolution_clock::now();
double dt = std::chrono::duration_cast<std::chrono::nanoseconds>(t2 - t1).count();
results.push_back( std::make_tuple( X.size(),
dt/1e6, is_fconvex ) );
if ( force_fconvexity && ! is_fconvex )
trace.
warning() <<
"Invalid computation of either FC* or full convexity !" << std::endl;
}
}
template <Dimension dim>
void
( std::vector< std::tuple< std::size_t, double, bool > >& results,
std::size_t nb_tries, std::size_t range )
{
typedef KhalimskySpaceND<dim,int64_t>
KSpace;
typedef PConvexity< Space > PConvexity;
DConvexity dconv( Point::diagonal( -1 ), Point::diagonal( range ) );
PConvexity pconv;
Domain domain( Point::diagonal( 0 ), Point::diagonal( range ) );
std::cout <<
"Computing " << nb_tries <<
" P-convexities in Z" <<
dim << std::endl;
for ( auto n = 0; n < nb_tries; ++n )
{
double filling_probability = 0.1 + 0.9 * double( n ) / double( nb_tries );
std::set< Point > S;
std::size_t nb_vertices
= std::size_t( filling_probability * ceil( pow( range,
dim ) ) );
for ( auto i = 0; i < nb_vertices; i++ ) {
for (
auto j = 0; j <
dim; j++ ) p[ j ] = rand() % range;
S.insert( p );
}
std::vector< Point > X( S.cbegin(), S.cend() );
std::chrono::high_resolution_clock::time_point
t1 = std::chrono::high_resolution_clock::now();
bool is_pconvex = pconv.isPConvex( X );
std::chrono::high_resolution_clock::time_point
t2 = std::chrono::high_resolution_clock::now();
double dt = std::chrono::duration_cast<std::chrono::nanoseconds>(t2 - t1).count();
results.push_back( std::make_tuple( X.size(),
dt/1e6, is_pconvex ) );
}
}
template <Dimension dim>
void
( std::vector< std::tuple< std::size_t, double, bool > >& results,
std::size_t nb_tries, std::size_t range )
{
typedef KhalimskySpaceND<dim,int64_t>
KSpace;
typedef PConvexity< Space > PConvexity;
DConvexity dconv( Point::diagonal( -1 ), Point::diagonal( range ) );
PConvexity pconv;
Domain domain( Point::diagonal( 0 ), Point::diagonal( range ) );
std::cout <<
"Computing " << nb_tries <<
" full convexities in Z" <<
dim << std::endl;
for ( auto n = 0; n < nb_tries; ++n )
{
double filling_probability = 0.1 + 0.9 * double( n ) / double( nb_tries );
std::set< Point > S;
std::size_t nb_vertices
= std::size_t( filling_probability * ceil( pow( range,
dim ) ) );
for ( auto i = 0; i < nb_vertices; i++ ) {
for (
auto j = 0; j <
dim; j++ ) p[ j ] = rand() % range;
S.insert( p );
}
std::vector< Point > X( S.cbegin(), S.cend() );
std::chrono::high_resolution_clock::time_point
t1 = std::chrono::high_resolution_clock::now();
std::chrono::high_resolution_clock::time_point
t2 = std::chrono::high_resolution_clock::now();
double dt = std::chrono::duration_cast<std::chrono::nanoseconds>(t2 - t1).count();
results.push_back( std::make_tuple( X.size(),
dt/1e6, is_fconvex ) );
}
}
template <Dimension dim>
void
( std::vector< std::tuple< std::size_t, double, bool > >& results,
std::size_t nb_tries, std::size_t range )
{
typedef KhalimskySpaceND<dim,int64_t>
KSpace;
typedef PConvexity< Space > PConvexity;
DConvexity dconv( Point::diagonal( -1 ), Point::diagonal( range ) );
PConvexity pconv;
Domain domain( Point::diagonal( 0 ), Point::diagonal( range ) );
std::cout <<
"Computing " << nb_tries <<
" full convexities (fast) in Z" <<
dim << std::endl;
for ( auto n = 0; n < nb_tries; ++n )
{
double filling_probability = 0.1 + 0.9 * double( n ) / double( nb_tries );
std::set< Point > S;
std::size_t nb_vertices
= std::size_t( filling_probability * ceil( pow( range,
dim ) ) );
for ( auto i = 0; i < nb_vertices; i++ ) {
for (
auto j = 0; j <
dim; j++ ) p[ j ] = rand() % range;
S.insert( p );
}
std::vector< Point > X( S.cbegin(), S.cend() );
std::chrono::high_resolution_clock::time_point
t1 = std::chrono::high_resolution_clock::now();
std::chrono::high_resolution_clock::time_point
t2 = std::chrono::high_resolution_clock::now();
double dt = std::chrono::duration_cast<std::chrono::nanoseconds>(t2 - t1).count();
results.push_back( std::make_tuple( X.size(),
dt/1e6, is_fconvex ) );
}
}
const std::vector< std::tuple< std::size_t, double, bool > >& results,
const std::string& fname )
{
std::ofstream output( fname );
output << "# Results of " << results.size() << " P-convexity computations in Z"
<< "# Card(X) time(ms) p-convex?" << std::endl;
for ( auto&& r : results )
output << std::get<0>( r ) << " " << std::get<1>( r ) << " " << std::get<2>( r )
<< std::endl;
output.close();
}
int main(
int argc,
char* argv[] )
{
if ( false )
{
std::vector< std::tuple< std::size_t, double, bool > > R2;
timingsPConvexity<2>( R2, 50, 3, 100, 0.5 );
timingsPConvexity<2>( R2, 50, 4, 200, 0.5 );
timingsPConvexity<2>( R2, 50, 5, 400, 0.5 );
timingsPConvexity<2>( R2, 50, 5, 600, 0.5 );
timingsPConvexity<2>( R2, 50, 5, 800, 0.5 );
timingsPConvexity<2>( R2, 25, 5,1200, 0.5 );
timingsPConvexity<2>( R2, 25, 5,2000, 0.5 );
}
if ( false )
{
std::vector< std::tuple< std::size_t, double, bool > > R3;
timingsPConvexity<3>( R3, 50, 3, 10, 0.5 );
timingsPConvexity<3>( R3, 50, 4, 20, 0.5 );
timingsPConvexity<3>( R3, 50, 5, 40, 0.5 );
timingsPConvexity<3>( R3, 50, 5, 80, 0.5 );
timingsPConvexity<3>( R3, 25, 5, 160, 0.5 );
timingsPConvexity<3>( R3, 25, 5, 320, 0.5 );
}
if ( false )
{
std::vector< std::tuple< std::size_t, double, bool > > R4;
timingsPConvexity<4>( R4, 50, 5, 10, 0.5 );
timingsPConvexity<4>( R4, 50, 5, 15, 0.5 );
timingsPConvexity<4>( R4, 50, 5, 20, 0.5 );
timingsPConvexity<4>( R4, 50, 5, 30, 0.5 );
timingsPConvexity<4>( R4, 25, 5, 40, 0.5 );
timingsPConvexity<4>( R4, 25, 5, 60, 0.5 );
timingsPConvexity<4>( R4, 15, 6, 80, 0.5 );
timingsPConvexity<4>( R4, 15, 6, 100, 0.5 );
timingsPConvexity<4>( R4, 15, 6, 120, 0.5 );
}
if ( false )
{
std::vector< std::tuple< std::size_t, double, bool > > R2;
timingsFullConvexity<2>( R2, 50, 3, 100, 0.5 );
timingsFullConvexity<2>( R2, 50, 4, 200, 0.5 );
timingsFullConvexity<2>( R2, 50, 5, 400, 0.5 );
timingsFullConvexity<2>( R2, 50, 5, 600, 0.5 );
timingsFullConvexity<2>( R2, 50, 5, 800, 0.5 );
timingsFullConvexity<2>( R2, 25, 5,1200, 0.5 );
timingsFullConvexity<2>( R2, 25, 5,2000, 0.5 );
}
if ( false )
{
std::vector< std::tuple< std::size_t, double, bool > > R3;
timingsFullConvexity<3>( R3, 50, 3, 10, 0.5 );
timingsFullConvexity<3>( R3, 50, 4, 20, 0.5 );
timingsFullConvexity<3>( R3, 50, 5, 40, 0.5 );
timingsFullConvexity<3>( R3, 50, 5, 80, 0.5 );
timingsFullConvexity<3>( R3, 25, 5, 160, 0.5 );
timingsFullConvexity<3>( R3, 25, 5, 320, 0.5 );
}
if ( false )
{
std::vector< std::tuple< std::size_t, double, bool > > R4;
timingsFullConvexity<4>( R4, 50, 5, 10, 0.5 );
timingsFullConvexity<4>( R4, 50, 5, 15, 0.5 );
timingsFullConvexity<4>( R4, 50, 5, 20, 0.5 );
timingsFullConvexity<4>( R4, 50, 5, 30, 0.5 );
timingsFullConvexity<4>( R4, 25, 5, 40, 0.5 );
timingsFullConvexity<4>( R4, 25, 5, 60, 0.5 );
timingsFullConvexity<4>( R4, 15, 6, 80, 0.5 );
timingsFullConvexity<4>( R4, 10, 6, 100, 0.5 );
timingsFullConvexity<4>( R4, 5, 6, 120, 0.5 );
}
if ( false )
{
std::vector< std::tuple< std::size_t, double, bool > > R2;
timingsFullConvexityFast<2>( R2, 50, 3, 100, 0.5 );
timingsFullConvexityFast<2>( R2, 50, 4, 200, 0.5 );
timingsFullConvexityFast<2>( R2, 50, 5, 400, 0.5 );
timingsFullConvexityFast<2>( R2, 50, 5, 600, 0.5 );
timingsFullConvexityFast<2>( R2, 50, 5, 800, 0.5 );
timingsFullConvexityFast<2>( R2, 25, 5,1200, 0.5 );
timingsFullConvexityFast<2>( R2, 25, 5,2000, 0.5 );
}
if ( false )
{
std::vector< std::tuple< std::size_t, double, bool > > R3;
timingsFullConvexityFast<3>( R3, 50, 3, 10, 0.5 );
timingsFullConvexityFast<3>( R3, 50, 4, 20, 0.5 );
timingsFullConvexityFast<3>( R3, 50, 5, 40, 0.5 );
timingsFullConvexityFast<3>( R3, 50, 5, 80, 0.5 );
timingsFullConvexityFast<3>( R3, 25, 5, 160, 0.5 );
timingsFullConvexityFast<3>( R3, 25, 5, 320, 0.5 );
}
if ( false )
{
std::vector< std::tuple< std::size_t, double, bool > > R4;
timingsFullConvexityFast<4>( R4, 50, 5, 10, 0.5 );
timingsFullConvexityFast<4>( R4, 50, 5, 15, 0.5 );
timingsFullConvexityFast<4>( R4, 50, 5, 20, 0.5 );
timingsFullConvexityFast<4>( R4, 50, 5, 30, 0.5 );
timingsFullConvexityFast<4>( R4, 25, 5, 40, 0.5 );
timingsFullConvexityFast<4>( R4, 25, 5, 60, 0.5 );
timingsFullConvexityFast<4>( R4, 15, 6, 80, 0.5 );
timingsFullConvexityFast<4>( R4, 10, 6, 100, 0.5 );
timingsFullConvexityFast<4>( R4, 5, 6, 120, 0.5 );
}
if ( false )
{
std::vector< std::tuple< std::size_t, double, bool > > R2;
timingsPConvexityNonConvex<2>( R2, 50, 100 );
timingsPConvexityNonConvex<2>( R2, 50, 200 );
timingsPConvexityNonConvex<2>( R2, 50, 400 );
timingsPConvexityNonConvex<2>( R2, 50, 600 );
timingsPConvexityNonConvex<2>( R2, 50, 800 );
timingsPConvexityNonConvex<2>( R2, 50, 1200 );
timingsPConvexityNonConvex<2>( R2, 50, 2000 );
}
if ( false )
{
std::vector< std::tuple< std::size_t, double, bool > > R3;
timingsPConvexityNonConvex<3>( R3, 50, 20 );
timingsPConvexityNonConvex<3>( R3, 50, 40 );
timingsPConvexityNonConvex<3>( R3, 50, 80 );
timingsPConvexityNonConvex<3>( R3, 50, 160 );
timingsPConvexityNonConvex<3>( R3, 50, 320 );
}
if ( false )
{
std::vector< std::tuple< std::size_t, double, bool > > R4;
timingsPConvexityNonConvex<4>( R4, 50, 10 );
timingsPConvexityNonConvex<4>( R4, 50, 20 );
timingsPConvexityNonConvex<4>( R4, 50, 30 );
timingsPConvexityNonConvex<4>( R4, 40, 40 );
timingsPConvexityNonConvex<4>( R4, 20, 50 );
}
if ( false )
{
std::vector< std::tuple< std::size_t, double, bool > > R2;
timingsFullConvexityNonConvex<2>( R2, 50, 100 );
timingsFullConvexityNonConvex<2>( R2, 50, 200 );
timingsFullConvexityNonConvex<2>( R2, 50, 400 );
timingsFullConvexityNonConvex<2>( R2, 50, 600 );
timingsFullConvexityNonConvex<2>( R2, 50, 800 );
timingsFullConvexityNonConvex<2>( R2, 50, 1200 );
timingsFullConvexityNonConvex<2>( R2, 50, 2000 );
}
if ( false )
{
std::vector< std::tuple< std::size_t, double, bool > > R3;
timingsFullConvexityNonConvex<3>( R3, 50, 20 );
timingsFullConvexityNonConvex<3>( R3, 50, 40 );
timingsFullConvexityNonConvex<3>( R3, 50, 80 );
timingsFullConvexityNonConvex<3>( R3, 40, 160 );
timingsFullConvexityNonConvex<3>( R3, 25, 320 );
}
if ( false )
{
std::vector< std::tuple< std::size_t, double, bool > > R4;
timingsFullConvexityNonConvex<4>( R4, 50, 10 );
timingsFullConvexityNonConvex<4>( R4, 50, 20 );
timingsFullConvexityNonConvex<4>( R4, 50, 30 );
timingsFullConvexityNonConvex<4>( R4, 40, 40 );
timingsFullConvexityNonConvex<4>( R4, 20, 50 );
}
if ( false )
{
std::vector< std::tuple< std::size_t, double, bool > > R2;
timingsFullConvexityFastNonConvex<2>( R2, 50, 100 );
timingsFullConvexityFastNonConvex<2>( R2, 50, 200 );
timingsFullConvexityFastNonConvex<2>( R2, 50, 400 );
timingsFullConvexityFastNonConvex<2>( R2, 50, 600 );
timingsFullConvexityFastNonConvex<2>( R2, 50, 800 );
timingsFullConvexityFastNonConvex<2>( R2, 50, 1200 );
timingsFullConvexityFastNonConvex<2>( R2, 50, 2000 );
}
if ( false )
{
std::vector< std::tuple< std::size_t, double, bool > > R3;
timingsFullConvexityFastNonConvex<3>( R3, 50, 20 );
timingsFullConvexityFastNonConvex<3>( R3, 50, 40 );
timingsFullConvexityFastNonConvex<3>( R3, 50, 80 );
timingsFullConvexityFastNonConvex<3>( R3, 40, 160 );
timingsFullConvexityFastNonConvex<3>( R3, 25, 320 );
}
if ( false )
{
std::vector< std::tuple< std::size_t, double, bool > > R4;
timingsFullConvexityFastNonConvex<4>( R4, 50, 10 );
timingsFullConvexityFastNonConvex<4>( R4, 50, 20 );
timingsFullConvexityFastNonConvex<4>( R4, 50, 30 );
timingsFullConvexityFastNonConvex<4>( R4, 40, 40 );
timingsFullConvexityFastNonConvex<4>( R4, 20, 50 );
}
return 0;
}
void getPoints(std::vector< Point > &pts) const
bool isFullyConvexFast(const PointRange &X) const
LatticePolytope CvxH(const PointRange &X) const
bool isFullyConvex(const PointRange &X, bool convex0=false) const
PointRange envelope(const PointRange &Z, EnvelopeAlgorithm algo=EnvelopeAlgorithm::DIRECT) const
double rand01()
[QuickHull3D-Includes]
DGtal is the top-level namespace which contains all DGtal functions and types.
void timingsFullConvexityNonConvex(std::vector< std::tuple< std::size_t, double, bool > > &results, std::size_t nb_tries, std::size_t range)
void timingsFullConvexity(std::vector< std::tuple< std::size_t, double, bool > > &results, std::size_t nb_tries, std::size_t nb_vertices, std::size_t range, double fconvexity_probability=0.5)
void outputResults(Dimension dim, const std::vector< std::tuple< std::size_t, double, bool > > &results, const std::string &fname)
void timingsPConvexityNonConvex(std::vector< std::tuple< std::size_t, double, bool > > &results, std::size_t nb_tries, std::size_t range)
void timingsFullConvexityFast(std::vector< std::tuple< std::size_t, double, bool > > &results, std::size_t nb_tries, std::size_t nb_vertices, std::size_t range, double fconvexity_probability=0.5)
void timingsPConvexity(std::vector< std::tuple< std::size_t, double, bool > > &results, std::size_t nb_tries, std::size_t nb_vertices, std::size_t range, double pconvexity_probability=0.5)
void timingsFullConvexityFastNonConvex(std::vector< std::tuple< std::size_t, double, bool > > &results, std::size_t nb_tries, std::size_t range)
int main(int argc, char **argv)
HyperRectDomain< Space > Domain
unsigned int dim(const Vector &z)