DGtal  1.2.beta
geometry/tools/exampleRationalBallDelaunay3D.cpp

Computation of the Delaunay complex of a set of rational points in 3D by Quick Hull algorithm.

./examples/geometry/tools/exampleRationalBallDelaunay3D 100 10 0.25 100

outputs

Compute convex hull in higher dimension
assign ridges/faces to cell and conversely
takes care of vertex positions
[ConvexCellComplex<3> #C=456 #F=955 #V=99 hasFaceGeometry ]
exampleRationalBallDelaunay3D-100.jpg
Delaunay cell decomposition of 100 randomly chosen points in a 3D ball with radius 10, with precision 100 and retraction 0.25
See also
QuickHull algorithm in arbitrary dimension for convex hull and Delaunay cell complex computation
#include "DGtal/base/Common.h"
#include "DGtal/kernel/PointVector.h"
#include "DGtal/shapes/SurfaceMesh.h"
#include "DGtal/io/writers/SurfaceMeshWriter.h"
#include "DGtal/geometry/volumes/ConvexityHelper.h"
using namespace DGtal;
using namespace DGtal::Z3i;
int main( int argc, char* argv[] )
{
int nb = argc > 1 ? atoi( argv[ 1 ] ) : 100; // nb points
double dR = argc > 2 ? atof( argv[ 2 ] ) : 10.; // radius of balla
double eps = argc > 3 ? atof( argv[ 3 ] ) : 0.1; // retraction
double precision = argc > 4 ? atof( argv[ 4 ] ) : 100.0; // precision
// (1) create range of random points in ball
std::vector< RealPoint > V;
const auto R2 = dR * dR;
for ( int i = 0; i < nb; ) {
RealPoint p( ( rand() / (double) RAND_MAX * 2.0 - 1.0 ) * dR,
( rand() / (double) RAND_MAX * 2.0 - 1.0 ) * dR,
( rand() / (double) RAND_MAX * 2.0 - 1.0 ) * dR );
if ( p.squaredNorm() <= R2 ) { V.push_back( p ); i++; }
}
// (2) compute convex hull
bool ok =
precision, true );
if ( ! ok )
{
trace.error() << "Input set of points is not full dimensional." << std::endl;
return 1;
}
dcomplex.requireFaceGeometry();
std::cout << dcomplex << std::endl;
// (3) build the mesh that is made of the exploded 3d cells
std::vector< RealPoint > positions;
std::vector< std::vector< Index > > facets;
Index idxv = 0;
for ( auto c = 0; c < dcomplex.nbCells(); ++c )
{
RealPoint b = dcomplex.cellBarycenter( c );
auto c_vtcs = dcomplex.cellVertices( c );
std::map< Index, Index > v2v;
for ( auto v : c_vtcs ) {
RealPoint x = dcomplex.position( v );
v2v[ v ] = idxv++;
positions.push_back( b + ( x - b ) * ( 1.0 - eps ) );
}
for ( const auto& f : dcomplex.cellFaces( c ) ) {
auto f_vtcs = dcomplex.faceVertices( f );
for ( auto& vertex : f_vtcs )
vertex = v2v[ vertex ];
facets.push_back( f_vtcs );
}
}
SMesh mesh( positions.cbegin(), positions.cend(),
facets.cbegin(), facets.cend() );
// (4) output result as OBJ file
std::ofstream out( "delaunay3d.obj" );
::writeOBJ( out, mesh );
out.close();
return 0;
}