Created
October 5, 2018 13:38
-
-
Save VincentRouvreau/54c5e0c2d1f9a6c94b0f15db37319775 to your computer and use it in GitHub Desktop.
gist_for_cgal_issue_3346
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h> | |
#include <CGAL/Periodic_3_Delaunay_triangulation_traits_3.h> | |
#include <CGAL/Periodic_3_Delaunay_triangulation_3.h> | |
#include <CGAL/Alpha_shape_3.h> | |
#include <CGAL/Alpha_shape_cell_base_3.h> | |
#include <CGAL/Alpha_shape_vertex_base_3.h> | |
#include <CGAL/iterator.h> | |
#include <CGAL/Random.h> | |
#include <CGAL/point_generators_3.h> | |
#include <fstream> | |
#include <cmath> | |
#include <string> | |
#include <tuple> | |
#include <map> | |
#include <utility> | |
#include <vector> | |
#include <cstdlib> | |
// Traits | |
using K = CGAL::Exact_predicates_inexact_constructions_kernel; | |
using PK = CGAL::Periodic_3_Delaunay_triangulation_traits_3<K>; | |
// Vertex type | |
using DsVb = CGAL::Periodic_3_triangulation_ds_vertex_base_3<>; | |
using Vb = CGAL::Triangulation_vertex_base_3<PK, DsVb>; | |
using AsVb = CGAL::Alpha_shape_vertex_base_3<PK, Vb>; | |
// Cell type | |
using DsCb = CGAL::Periodic_3_triangulation_ds_cell_base_3<>; | |
using Cb = CGAL::Triangulation_cell_base_3<PK, DsCb>; | |
using AsCb = CGAL::Alpha_shape_cell_base_3<PK, Cb>; | |
using Tds = CGAL::Triangulation_data_structure_3<AsVb, AsCb>; | |
using P3DT3 = CGAL::Periodic_3_Delaunay_triangulation_3<PK, Tds>; | |
using Alpha_shape_3 = CGAL::Alpha_shape_3<P3DT3>; | |
using Point = PK::Point_3; | |
// filtration with alpha values needed type definition | |
using Alpha_value_type = Alpha_shape_3::FT; | |
using Object = CGAL::Object; | |
using Dispatch = | |
CGAL::Dispatch_output_iterator<CGAL::cpp11::tuple<Object, Alpha_value_type>, | |
CGAL::cpp11::tuple<std::back_insert_iterator<std::vector<Object> >, | |
std::back_insert_iterator<std::vector<Alpha_value_type> > > >; | |
using Cell_handle = Alpha_shape_3::Cell_handle; | |
using Facet = Alpha_shape_3::Facet; | |
using Edge_3 = Alpha_shape_3::Edge; | |
using Vertex_handle = Alpha_shape_3::Vertex_handle; | |
using Vertex_list = std::vector<Alpha_shape_3::Vertex_handle>; | |
int main(int argc, char **argv) { | |
typedef CGAL::Creator_uniform_3<double, Point> Creator; | |
CGAL::Random random(7); | |
CGAL::Random_points_in_cube_3<Point, Creator> in_cube(1, random); | |
std::vector<Point> pts; | |
// Generating 50 random points | |
for (int i=0 ; i < 50 ; i++) { | |
Point p = *in_cube++; | |
pts.push_back(p); | |
} | |
for (int i=0 ; i < 5 ; i++) { | |
std::ofstream out(std::to_string(i) + std::string(".txt")); | |
// I maintain a map from Alpha_shape_3::Vertex_handle to my own vertex handles | |
std::map<Alpha_shape_3::Vertex_handle, int> map_cgal_simplex_tree; | |
// Define the periodic cube | |
P3DT3 pdt(PK::Iso_cuboid_3(-1,-1,-1,1,1,1)); | |
// Heuristic for inserting large point sets (if pts is reasonably large) | |
pdt.insert(pts.begin(), pts.end(), true); | |
// As pdt won't be modified anymore switch to 1-sheeted cover if possible | |
if (pdt.is_triangulation_in_1_sheet()) { | |
pdt.convert_to_1_sheeted_covering(); | |
} else { | |
std::cerr << "ERROR: we were not able to construct a triangulation within a single periodic domain." << std::endl; | |
exit(-1); | |
} | |
std::cout << "============================ Periodic Delaunay computed." << std::endl; | |
// alpha shape construction from points. CGAL has a strange behavior in REGULARIZED mode. This is the default mode | |
// Maybe need to set it to GENERAL mode | |
Alpha_shape_3 as(pdt, 0, Alpha_shape_3::GENERAL); | |
// filtration with alpha values from alpha shape | |
std::vector<Object> the_objects; | |
std::vector<Alpha_value_type> the_alpha_values; | |
Dispatch disp = CGAL::dispatch_output<Object, Alpha_value_type>(std::back_inserter(the_objects), | |
std::back_inserter(the_alpha_values)); | |
as.filtration_with_alpha_values(disp); | |
// Loop on objects vector | |
for (auto object_iterator : the_objects) { | |
Vertex_list vertex_list; | |
// Retrieve Alpha shape vertex list from object | |
if (const Cell_handle *ch = CGAL::object_cast<Cell_handle>(&object_iterator)) { | |
for (auto i = 0; i < 4; i++) { | |
vertex_list.push_back((*ch)->vertex(i)); | |
} | |
} else if (const Facet *facet = CGAL::object_cast<Facet>(&object_iterator)) { | |
for (auto i = 0; i < 4; i++) { | |
if (facet->second != i) { | |
vertex_list.push_back(facet->first->vertex(i)); | |
} | |
} | |
} else if (const Edge_3 *edge = CGAL::object_cast<Edge_3>(&object_iterator)) { | |
for (auto i : {edge->second, edge->third}) { | |
vertex_list.push_back(edge->first->vertex(i)); | |
} | |
} else if (const Vertex_handle *vh = CGAL::object_cast<Vertex_handle>(&object_iterator)) { | |
vertex_list.push_back(*vh); | |
} | |
std::vector<int> the_simplex; | |
// Construction of the vector of simplex_tree vertex from list of alpha_shapes vertex | |
for (auto the_alpha_shape_vertex : vertex_list) { | |
std::map<Alpha_shape_3::Vertex_handle, int>::iterator the_map_iterator = map_cgal_simplex_tree.find(the_alpha_shape_vertex); | |
if (the_map_iterator == map_cgal_simplex_tree.end()) { | |
// alpha shape not found | |
int vertex = map_cgal_simplex_tree.size(); | |
#ifdef DEBUG_TRACES | |
//std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] not found - insert " << vertex << std::endl; | |
#endif // DEBUG_TRACES | |
map_cgal_simplex_tree.emplace(the_alpha_shape_vertex, vertex); | |
the_simplex.push_back(vertex); | |
} else { | |
// alpha shape found | |
int vertex = the_map_iterator->second; | |
#ifdef DEBUG_TRACES | |
//std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] found in " << vertex << std::endl; | |
#endif // DEBUG_TRACES | |
the_simplex.push_back(vertex); | |
} | |
} | |
for (auto vertex : the_simplex) | |
out << vertex << ","; | |
out << std::endl; | |
} | |
out.close(); | |
} | |
return 0; | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment