DGtal  1.4.beta
checkLatticeBallQuickHull.cpp
Go to the documentation of this file.
1 
61 #include <cstdlib>
62 #include <iostream>
63 #include <chrono>
64 #include "DGtal/base/Common.h"
65 #include "DGtal/geometry/tools/QuickHull.h"
66 
67 template <typename Point>
68 std::vector< Point >
69 randomPointsInBall( int nb, double R )
70 {
71  std::vector< Point > V;
72  DGtal::int64_t iR = DGtal::int64_t( round( R ) );
73  Point c = Point::diagonal( iR );
74  double R2 = (double) R * (double) R;
75  for ( int i = 0; i < nb; ) {
76  Point p;
77  for ( DGtal::Dimension k = 0; k < Point::dimension; ++k )
78  p[ k ] = DGtal::int64_t( round( (double) rand() * 2.0 * R / (double) RAND_MAX ));
79  if ( ( p - c ).squaredNorm() < R2 ) { V.push_back( p - c ); i++; }
80  }
81  return V;
82 }
83 
84 template <typename Point>
85 std::vector< Point >
86 randomPointsInBallBigInteger( int nb, double R )
87 {
88  std::vector< Point > V;
89  DGtal::int64_t iR = DGtal::int64_t( round( R ) );
91  Point c = Point::diagonal( Converter::cast( iR ) );
92  double R2 = (double) R * (double) R;
93  for ( int i = 0; i < nb; ) {
94  Point p;
95  for ( DGtal::Dimension k = 0; k < Point::dimension; ++k )
96  p[ k ] = Converter::cast( DGtal::int64_t( round( (double) rand() * 2.0 * R / (double) RAND_MAX )) );
97  if ( ( p - c ).squaredNorm() < R2 ) { V.push_back( p - c ); i++; }
98  }
99  return V;
100 }
101 
102 template < DGtal::Dimension dim, typename Integer >
103 bool
104 checkQuickHull( int nb, double R )
105 {
107  typedef DGtal::QuickHull< Kernel > ConvexHull;
108  typedef typename ConvexHull::Point Point;
109 
110  const auto V = randomPointsInBall< Point >( nb, R );
111  ConvexHull hull;
112  hull.setInput( V );
113  hull.computeConvexHull();
114  std::cout << "#points=" << hull.nbPoints()
115  << " #vertices=" << hull.nbVertices()
116  << " #facets=" << hull.nbFacets() << std::endl;
117  double total_time = 0;
118  std::for_each( hull.timings.cbegin(), hull.timings.cend(),
119  [&total_time] ( double t ) { total_time += t; } );
120  std::cout << "purge duplicates= " << round(hull.timings[ 0 ]) << " ms." << std::endl;
121  std::cout << "init simplex = " << round(hull.timings[ 1 ]) << " ms." << std::endl;
122  std::cout << "quickhull core = " << round(hull.timings[ 2 ]) << " ms." << std::endl;
123  std::cout << "compute vertices= " << round(hull.timings[ 3 ]) << " ms." << std::endl;
124  std::cout << "total time = " << round(total_time) << " ms." << std::endl;
125  std::cout << "Checking hull ... " << std::endl;
126  auto start = std::chrono::steady_clock::now();
127  bool ok = hull.check();
128  auto end = std::chrono::steady_clock::now();
129  std::chrono::duration<double> elapsed_seconds = end-start;
130  std::cout << " ... in " << elapsed_seconds.count() << "s"
131  << " => " << ( ok ? "OK" : "ERROR" ) << std::endl;
132  return ok;
133 }
134 
135 template < DGtal::Dimension dim >
136 bool
137 checkQuickHullBigInteger( int nb, double R )
138 {
140  typedef DGtal::QuickHull< Kernel > ConvexHull;
141  typedef typename ConvexHull::Point Point;
142 
143  const auto V = randomPointsInBallBigInteger< Point >( nb, R );
144  ConvexHull hull;
145  hull.setInput( V );
146  hull.computeConvexHull();
147  std::cout << "#points=" << hull.nbPoints()
148  << " #vertices=" << hull.nbVertices()
149  << " #facets=" << hull.nbFacets() << std::endl;
150  double total_time = 0;
151  std::for_each( hull.timings.cbegin(), hull.timings.cend(),
152  [&total_time] ( double t ) { total_time += t; } );
153  std::cout << "purge duplicates= " << round(hull.timings[ 0 ]) << " ms." << std::endl;
154  std::cout << "init simplex = " << round(hull.timings[ 1 ]) << " ms." << std::endl;
155  std::cout << "quickhull core = " << round(hull.timings[ 2 ]) << " ms." << std::endl;
156  std::cout << "compute vertices= " << round(hull.timings[ 3 ]) << " ms." << std::endl;
157  std::cout << "total time = " << round(total_time) << " ms." << std::endl;
158  std::cout << "Checking hull ... " << std::endl;
159  auto start = std::chrono::steady_clock::now();
160  bool ok = hull.check();
161  auto end = std::chrono::steady_clock::now();
162  std::chrono::duration<double> elapsed_seconds = end-start;
163  std::cout << " ... in " << elapsed_seconds.count() << "s"
164  << " => " << ( ok ? "OK" : "ERROR" ) << std::endl;
165  return ok;
166 }
167 
168 int main( int argc, char* argv[] )
169 {
170  int dim = argc > 1 ? atoi( argv[ 1 ] ) : 3; // dimension
171  int nb = argc > 2 ? atoi( argv[ 2 ] ) : 1000; // nb points
172  double R = argc > 3 ? atof( argv[ 3 ] ) : 100.0; // radius of ball
173  std::string i = argc > 4 ? argv[ 4 ] : "int64"; // type for internal integers
174  bool ok = true;
175  if ( ( i != "int64" ) && ( i != "bigint" ) && ( i != "allbigint" ) )
176  {
177  DGtal::trace.error() << "Integer type in {int64,bigint,allbigint}" << std::endl;
178  ok = false;
179  }
180  if ( ( dim < 2 ) || ( dim > 6 ) )
181  {
182  DGtal::trace.error() << "Dimension must be in {2,3,4,5,6}" << std::endl;
183  ok = false;
184  }
185  if ( ! ok ) return 1;
186  if ( i == "bigint" )
187  {
188  switch( dim ) {
189  case 2 : ok = checkQuickHull< 2, DGtal::BigInteger >( nb, R ); break;
190  case 3 : ok = checkQuickHull< 3, DGtal::BigInteger >( nb, R ); break;
191  case 4 : ok = checkQuickHull< 4, DGtal::BigInteger >( nb, R ); break;
192  case 5 : ok = checkQuickHull< 5, DGtal::BigInteger >( nb, R ); break;
193  case 6 : ok = checkQuickHull< 6, DGtal::BigInteger >( nb, R ); break;
194  }
195  }
196  else if ( i == "int64" )
197  {
198  switch( dim ) {
199  case 2 : ok = checkQuickHull< 2, DGtal::int64_t >( nb, R ); break;
200  case 3 : ok = checkQuickHull< 3, DGtal::int64_t >( nb, R ); break;
201  case 4 : ok = checkQuickHull< 4, DGtal::int64_t >( nb, R ); break;
202  case 5 : ok = checkQuickHull< 5, DGtal::int64_t >( nb, R ); break;
203  case 6 : ok = checkQuickHull< 6, DGtal::int64_t >( nb, R ); break;
204  }
205  }
206  else if ( i == "allbigint" )
207  {
208  switch( dim ) {
209  case 2 : ok = checkQuickHullBigInteger< 2 >( nb, R ); break;
210  case 3 : ok = checkQuickHullBigInteger< 3 >( nb, R ); break;
211  case 4 : ok = checkQuickHullBigInteger< 4 >( nb, R ); break;
212  case 5 : ok = checkQuickHullBigInteger< 5 >( nb, R ); break;
213  case 6 : ok = checkQuickHullBigInteger< 6 >( nb, R ); break;
214  }
215  }
216  return ok ? 0 : 1;
217 }
218 
bool checkQuickHull(int nb, double R)
int main(int argc, char *argv[])
bool checkQuickHullBigInteger(int nb, double R)
std::vector< Point > randomPointsInBall(int nb, double R)
std::vector< Point > randomPointsInBallBigInteger(int nb, double R)
std::ostream & error()
boost::int64_t int64_t
signed 94-bit integer.
Definition: BasicTypes.h:74
DGtal::uint32_t Dimension
Definition: Common.h:136
Trace trace
Definition: Common.h:153
Aim: a geometric kernel to compute the convex hull of digital points with integer-only arithmetic.
----------— INTEGER/POINT CONVERSION SERVICES -----------------—
Aim: Implements the quickhull algorithm by Barber et al. , a famous arbitrary dimensional convex hull...
Definition: QuickHull.h:140
MyPointD Point
Definition: testClone2.cpp:383
unsigned int dim(const Vector &z)