exampleGenericLatticeConvexHull3D.cpp

Go to the documentation of this file.

 
#include <iostream>
#include <vector>
#include <random>
#include <algorithm>
#include <polyscope/polyscope.h>
#include <polyscope/surface_mesh.h>
#include <polyscope/point_cloud.h>
#include <polyscope/curve_network.h>
 
#include "DGtal/base/Common.h"
#include "DGtal/helpers/StdDefs.h"
#include "DGtal/geometry/tools/GenericLatticeConvexHull.h"
 
 
using namespace DGtal;
 
 
//Polyscope global
polyscope::PointCloud   *psPoints;
polyscope::PointCloud   *psVertices;
polyscope::PointCloud   *psBoundary0;
polyscope::CurveNetwork *psBoundary1;
polyscope::SurfaceMesh  *psBoundary2;
 
std::random_device rd;
std::mt19937 g(rd());
 
template < typename Point >
static
std::vector< Point >
makeRandomLatticePointsFromDirVectors( Point A, const std::vector< Point>& V,
                                       int nb, double radius, int amplitude, int aff_dim )
{
  std::uniform_int_distribution<int> U(-amplitude, amplitude);
  std::vector< Point > P;
  int m = std::min( aff_dim, (int) V.size() );
  for ( auto k = 0; P.size() < nb && k < 100000; k++ )
    {
      Point B = A;
      for ( auto i = 0; i < m; i++ )
        {
          int l = U( g );
          B += l * V[ i ];
        }
      if ( (B-A).norm() <= radius )
        P.push_back( B );
    }
  std::shuffle( P.begin(), P.end(), g );
  return P;
}
 
int main( int argc, char* argv[] )
{
  typedef GenericLatticeConvexHull< 3, int > QHull;
  typedef SpaceND< 3, int >                  Space;
  typedef Space::Point                       Point;
 
  std::cout << "Usage: " << argv[ 0 ] << " [R=30] [N=30] [D=2]\n";
  std::cout << "Computes the convex hull of N points within a ball of radius R, these points belonging to a lattice of chosen dimension D.\n";
  double radius = argc > 1 ? atof( argv[ 1 ] ) : 30.0;
  int    nb     = argc > 2 ? atoi( argv[ 2 ] ) : 30;
  int    adim   = argc > 3 ? atoi( argv[ 3 ] ) : 2;
  if ( nb < 0 ) return 1;
  if ( adim < 0 || adim > 3 ) return 1;
 
  // Create points
  std::vector< Point > L = { Point{ 4, 1, -3 }, Point{ 0, 2, 5 }, Point{ -1, -3, 5 } };
  std::vector< Point > X
    = makeRandomLatticePointsFromDirVectors( Point(1,2,-1),
                                             L, nb,
                                             radius,
                                             int( round( radius+0.5 ) ),
                                             adim );
 
  // Compute convex hull
  QHull hull;
  bool ok = hull.compute( X );
  std:: cout << ( ok ? "[PASSED]" : "[FAILED]" ) << " hull=" << hull << "\n";
  // Initialize polyscope
  polyscope::init();
  psPoints   = polyscope::registerPointCloud( "Points",   X );
  psVertices = polyscope::registerPointCloud( "Vertices", hull.positions );
  if ( hull.affine_dimension <= 1 ) // 1 or 2 points
    {
      // no facets
      psBoundary0 = polyscope::registerPointCloud( "Convex hull bdy dim=0",
                                                   hull.positions );
    }
  else if ( hull.affine_dimension == 2 ) // 2D
    {
      // facets are edges (and implicitly converted by polyscope)
      psBoundary1 = polyscope::registerCurveNetwork( "Convex hull bdy dim=1",
                                                     hull.positions, hull.facets );
      psBoundary0 = polyscope::registerPointCloud( "Projected points",
                                                   hull.projected_points );
      std::set<Point> S( hull.projected_points.cbegin(),
                         hull.projected_points.cend() );
      std::cout << "Projection basis=[ " << hull.affine_basis.basis()[0]
                << "," << hull.affine_basis.basis()[1] << " ]"
                << " d=" << hull.projected_dilation << "\n";
    }
  else if ( hull.affine_dimension == 3 ) // 3D
    {
      // facets are polygons
      psBoundary2 = polyscope::registerSurfaceMesh("Convex hull bdy dim=2",
                                                   hull.positions, hull.facets );
    }
  std::cout << "      dilation=" << hull.projected_dilation
            << " => counting of lattice points is "
            << (hull.projected_dilation == 1 ? "correct" : "INCORRECT") << ".\n"; 
  std::cout << "     #(P ∩ Z3)=" << hull.count() << "\n";
  std::cout << "#(Int(P) ∩ Z3)=" << hull.countInterior() << "\n";
  std::cout << " #(Bd(P) ∩ Z3)=" << hull.countBoundary() << "\n";
  polyscope::show();
  return EXIT_SUCCESS;
  
}