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"
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 typename KSpace::Point
Point;
typedef typename KSpace::Space
Space;
DConvexity dconv( Point::diagonal( -1 ), Point::diagonal( range ) );
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();
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 typename KSpace::Point
Point;
typedef typename KSpace::Space
Space;
DConvexity dconv( Point::diagonal( -1 ), Point::diagonal( range ) );
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 typename KSpace::Point
Point;
typedef typename KSpace::Space
Space;
DConvexity dconv( Point::diagonal( -1 ), Point::diagonal( range ) );
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 typename KSpace::Point
Point;
typedef typename KSpace::Space
Space;
DConvexity dconv( Point::diagonal( -1 ), Point::diagonal( range ) );
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();
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 typename KSpace::Point
Point;
typedef typename KSpace::Space
Space;
DConvexity dconv( Point::diagonal( -1 ), Point::diagonal( range ) );
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 typename KSpace::Point
Point;
typedef typename KSpace::Space
Space;
DConvexity dconv( Point::diagonal( -1 ), Point::diagonal( range ) );
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;
}
if ( false )
{
std::vector< std::tuple< std::size_t, double, bool > > R3;
}
if ( false )
{
std::vector< std::tuple< std::size_t, double, bool > > R4;
}
if ( false )
{
std::vector< std::tuple< std::size_t, double, bool > > R2;
}
if ( false )
{
std::vector< std::tuple< std::size_t, double, bool > > R3;
}
if ( false )
{
std::vector< std::tuple< std::size_t, double, bool > > R4;
}
if ( false )
{
std::vector< std::tuple< std::size_t, double, bool > > R2;
}
if ( false )
{
std::vector< std::tuple< std::size_t, double, bool > > R3;
}
if ( false )
{
std::vector< std::tuple< std::size_t, double, bool > > R4;
}
if ( false )
{
std::vector< std::tuple< std::size_t, double, bool > > R2;
}
if ( false )
{
std::vector< std::tuple< std::size_t, double, bool > > R3;
}
if ( false )
{
std::vector< std::tuple< std::size_t, double, bool > > R4;
}
if ( false )
{
std::vector< std::tuple< std::size_t, double, bool > > R2;
}
if ( false )
{
std::vector< std::tuple< std::size_t, double, bool > > R3;
}
if ( false )
{
std::vector< std::tuple< std::size_t, double, bool > > R4;
}
if ( false )
{
std::vector< std::tuple< std::size_t, double, bool > > R2;
}
if ( false )
{
std::vector< std::tuple< std::size_t, double, bool > > R3;
}
if ( false )
{
std::vector< std::tuple< std::size_t, double, bool > > R4;
}
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
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...
Aim: A class to check if digital sets are P-convex. The P-convexity is defined as follows: A digital ...
bool isPConvex(const std::vector< Point > &X) const
double rand01()
[QuickHull3D-Includes]
DGtal is the top-level namespace which contains all DGtal functions and types.
DGtal::uint32_t Dimension
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)
HyperRectDomain< Space > Domain