geometry/tools/exampleLatticeBallDelaunay2D.cpp

Computation of the Delaunay complex of a set of lattice points in 2D by Quick Hull algorithm.

# 1000000 points in digital ball of radius 1000
./examples/geometry/tools/exampleLatticeBallDelaunay2D 1000000 1000

outputs

#points=999908 #vertices=4480 #facets=8947
purge duplicates= 281 ms.
init simplex    = 24 ms.
quickhull core  = 336 ms.
compute vertices= 24 ms.
total time      = 665 ms.

See also

QuickHull algorithm in arbitrary dimension for convex hull and Delaunay cell complex computation

 
#include "DGtal/base/Common.h"
#include "DGtal/geometry/tools/QuickHull.h"
#include "DGtal/shapes/SurfaceMesh.h"
#include "DGtal/io/writers/SurfaceMeshWriter.h"
 
using namespace DGtal::Z2i;
int main( int argc, char* argv[] )
{
  int    nb = argc > 1 ? atoi( argv[ 1 ] ) : 100;  // nb points
  double dR = argc > 2 ? atof( argv[ 2 ] ) : 10.0; // radius of ball
  // (0) typedefs
  typedef DGtal::DelaunayIntegralKernel< 2 > Kernel2D;
  typedef DGtal::QuickHull< Kernel2D >       Delaunay2D;
  // (1) create range of random points in ball
  std::vector< Point > V;
  const double R2 = dR * dR;
  const int    R  = (int) ceil( dR );
  for ( int i = 0; i < nb; ) {
    Point p( rand() % (2*R+1) - R, rand() % (2*R+1) - R, rand() % (2*R+1) - R );
    if ( p.squaredNorm() < R2 ) { V.push_back( p ); i++; }
  }
  // (2) compute convex hull
  Delaunay2D hull;
  hull.setInput( V );
  hull.computeConvexHull();
  std::cout << "#points="    << hull.nbPoints()
            << " #vertices=" << hull.nbVertices()
            << " #facets="   << hull.nbFacets() << std::endl;
  double total_time = 0;
  std::for_each( hull.timings.cbegin(), hull.timings.cend(),
                 [&total_time] ( double t ) { total_time += t; } );
  std::cout << "purge duplicates= " << round(hull.timings[ 0 ]) << " ms." << std::endl;
  std::cout << "init simplex    = " << round(hull.timings[ 1 ]) << " ms." << std::endl;
  std::cout << "quickhull core  = " << round(hull.timings[ 2 ]) << " ms." << std::endl;
  std::cout << "compute vertices= " << round(hull.timings[ 3 ]) << " ms." << std::endl;
  std::cout << "total time      = " << round(total_time) << " ms." << std::endl;
  // (3) build mesh
  std::vector< RealPoint > positions;
  hull.getVertexPositions( positions );
  std::vector< std::vector< std::size_t > > facets;
  hull.getFacetVertices( facets );
  // Finite facets precede infinite facets => keep only finite facets
  facets.resize( hull.nbFiniteFacets() ); 
  // To viusalize the result, we build a surface mesh in R3 lying in
  // the plane z=0 and composed of the Delaunay cells.
  typedef DGtal::SpaceND< 3, int > Z3;
  std::vector< Z3::RealPoint > positions3d;
  for ( auto p : positions )
    positions3d.push_back( Z3::RealPoint( p[ 0 ], p[ 1 ], 0.0  ) );
  typedef DGtal::SurfaceMesh< Z3::RealPoint, Z3::RealVector> SMesh;
  SMesh mesh( positions3d.cbegin(), positions3d.cend(),
              facets.cbegin(), facets.cend() );
  // (4) output result as OBJ file
  std::ofstream out( "delaunay.obj" );
  DGtal::SurfaceMeshWriter< Z3::RealPoint, Z3::RealVector >::writeOBJ( out, mesh );
  out.close();
  return 0;
}