delaunay_4.cpp

#include <CGAL/Epick_d.h>
#include <CGAL/point_generators_d.h>
#include <CGAL/Delaunay_triangulation.h>
#include <CGAL/algorithm.h>
#include <CGAL/Timer.h>
#include <CGAL/assertions.h>
#include <CGAL/property_map.h>
#include <CGAL/Spatial_sort_traits_adapter_d.h>

#include <iostream>
#include <iterator>
#include <vector>



template<int dim>
void compute_triangulation(std::vector<std::vector<double>> points_) {
	typedef CGAL::Epick_d< CGAL::Dimension_tag<dim> >  K;
	typedef CGAL::Triangulation_data_structure<typename K::Dimension, CGAL::Triangulation_vertex<K,int>, CGAL::Triangulation_full_cell<K,int>> TDS;
	typedef CGAL::Delaunay_triangulation<K,TDS> T;
	std::vector<typename T::Point> points;
	points.reserve(points_.size());
	for (int i = 0; i < points_.size(); i++) {
		points.push_back(T::Point(points_[i].begin(), points_[i].end()));
	}

	T t(dim);
	CGAL::Timer cost;  
	cost.start();
	std::cout << "Delaunay triangulation of " << points.size() << " points in dim " << dim << std::endl << std::flush;

	std::vector<std::size_t> indices;
	std::vector<int> infos;
	for (int i = 0; i < points.size(); i++) {
		infos.push_back(i);
		indices.push_back(i);
	}

	typedef typename CGAL::Pointer_property_map<typename T::Point>::type Pmap;
	typedef CGAL::Spatial_sort_traits_adapter_d<K, Pmap> Search_traits;

	spatial_sort(indices.begin(), indices.end(),
		Search_traits(make_property_map(points), t.geom_traits()));

	typename T::Full_cell_handle hint;
	for (int i = 0; i < points.size(); i++) {
		auto vertex=t.insert(points[indices[i]], hint);
		vertex->data() = indices[i];
		hint = vertex->full_cell();
	}


	std::cout << "Done in " << cost.time() << " seconds." << std::endl;
	std::cout << "Number of simplices: " << t.number_of_finite_full_cells() << std::endl;


}


int main(int argc, char **argv)
{
	int n = 20;
	if( argc > 1 ) n = atoi(argv[1]); // number of points

	std::vector<std::vector<double>> points;
	points.reserve(n*n*n);
	for (int i=0;i<n;i++){
		for (int j=0;j<n;j++){
			for (int k=0;k<n;k++){
				for (int l=0;l<n;l++) {
					points.push_back(std::vector<double>{(double)i, (double)j, (double)k, (double)l});
				}
			}
		}
	}

	compute_triangulation<4>(points);
  
	return 0;
}