16#include <gudhi/Clock.h>
17#include <gudhi/Skeleton_blocker/Skeleton_blocker_simple_geometric_traits.h>
18#include <gudhi/Skeleton_blocker_geometric_complex.h>
19#include <gudhi/Off_reader.h>
21#include <CGAL/Euclidean_distance.h>
28#include "utils/UI_utils.h"
29#include "utils/Lloyd_builder.h"
30#include "utils/Rips_builder.h"
31#include "utils/K_nearest_builder.h"
32#include "utils/Vertex_collapsor.h"
33#include "utils/Edge_collapsor.h"
34#include "utils/Edge_contractor.h"
35#include "utils/Persistence_compute.h"
36#include "utils/Critical_points.h"
37#include "utils/Is_manifold.h"
39#include "Complex_typedefs.h"
41template<
typename Complex>
42class CGAL_geometric_flag_complex_wrapper {
45 typedef typename Complex::Point Point;
47 const bool load_only_points_;
50 CGAL_geometric_flag_complex_wrapper(
Complex& complex,
bool load_only_points =
false) :
52 load_only_points_(load_only_points) { }
54 void init(
int dim,
int num_vertices,
int num_max_faces,
int num_edges)
const { }
56 void point(
const std::vector<double>& coords) {
57 Point p(coords.size(), coords.begin(), coords.end());
61 void maximal_face(std::vector<int> vertices) {
62 if (!load_only_points_) {
64 for (std::size_t i = 0; i < vertices.size(); ++i)
65 for (std::size_t j = i + 1; j < vertices.size(); ++j)
78 Model() : complex_() { }
81 void off_file_open(
const std::string& name_file) {
82 UIDBGMSG(
"load off file", name_file);
84 CGAL_geometric_flag_complex_wrapper<Complex> read_wraper(complex_,
false);
85 Gudhi::read_off(name_file, read_wraper);
88 void off_points_open(
const std::string& name_file) {
89 UIDBGMSG(
"load off points", name_file);
91 CGAL_geometric_flag_complex_wrapper<Complex> read_wraper(complex_,
true);
92 Gudhi::read_off(name_file, read_wraper);
95 void off_file_save(
const std::string& name_file) {
96 UIDBG(
"save off file");
97 UIDBG(
"save off off_points_save");
98 std::ofstream file(name_file);
101 file << complex_.num_vertices() <<
" " << complex_.num_edges() <<
" 0\n";
103 const auto& pt(complex_.
point(v));
104 for (
auto it = pt.cartesian_begin(); it != pt.cartesian_end(); ++it)
112 std::cerr <<
"Could not open file " << name_file << std::endl;
116 void off_points_save(
const std::string& name_file) {
117 UIDBG(
"save off off_points_save");
118 std::ofstream file(name_file);
119 if (file.is_open()) {
121 file << complex_.num_vertices() <<
" 0 0\n";
123 const auto& pt(complex_.
point(v));
124 for (
auto it = pt.cartesian_begin(); it != pt.cartesian_end(); ++it)
130 std::cerr <<
"Could not open file " << name_file << std::endl;
135 void uniform_noise(
double amplitude) {
138 complex_.
point(v) = add_uniform_noise(complex_.
point(v), amplitude);
142 Point add_uniform_noise(
const Point& point,
double amplitude) {
143 std::vector<double> new_point(point.dimension());
144 for (
int i = 0; i < point.dimension(); ++i) {
145 new_point[i] = point[i] + (rand() % 2 - .5) * amplitude;
147 return Point(point.dimension(), new_point.begin(), new_point.end());
151 void lloyd(
int num_iterations,
int num_closest_neighbors) {
156 double squared_eucl_distance(
const Point& p1,
const Point& p2)
const {
157 return Geometry_trait::Squared_distance_d()(p1, p2);
162 void build_rips(
double alpha) {
164 Rips_builder<Complex> rips_builder(complex_, alpha);
167 void build_k_nearest_neighbors(
unsigned k) {
168 UIDBG(
"build_k_nearest");
170 K_nearest_builder<Complex> k_nearest_builder(complex_, k);
173 void build_delaunay() {
174 UIDBG(
"build_delaunay");
178 void contract_edges(
unsigned num_contractions) {
181 std::clog <<
"Time to simplify: " << c.num_seconds() <<
"s" << std::endl;
184 void collapse_vertices(
unsigned num_collapses) {
185 auto old_num_vertices = complex_.num_vertices();
187 UIDBGMSG(
"num vertices collapsed:", old_num_vertices - complex_.num_vertices());
190 void collapse_edges(
unsigned num_collapses) {
194 void show_graph_stats() {
195 std::clog <<
"++++++ Graph stats +++++++" << std::endl;
196 std::clog <<
"Num vertices : " << complex_.num_vertices() << std::endl;
197 std::clog <<
"Num edges : " << complex_.num_edges() << std::endl;
199 std::clog <<
"Min/avg/max degree : " << min_degree() <<
"/" << avg_degree() <<
"/" << max_degree() << std::endl;
202 std::clog <<
"+++++++++++++++++++++++++" << std::endl;
206 int min_degree()
const {
207 int res = (std::numeric_limits<int>::max)();
209 res = (std::min)(res, complex_.
degree(v));
213 int max_degree()
const {
216 res = (std::max)(res, complex_.
degree(v));
220 int avg_degree()
const {
223 res += complex_.
degree(v);
224 return res / complex_.num_vertices();
228 void show_complex_stats() {
229 std::clog <<
"++++++ Mesh stats +++++++" << std::endl;
230 std::clog <<
"Num vertices : " << complex_.num_vertices() << std::endl;
231 std::clog <<
"Num edges : " << complex_.num_edges() << std::endl;
233 std::clog <<
"+++++++++++++++++++++++++" << std::endl;
236 void show_complex_dimension() {
237 unsigned num_simplices = 0;
243 dimension = (std::max)(s.dimension(), dimension);
244 if (s.dimension() % 2 == 0)
250 std::clog <<
"++++++ Mesh dimension +++++++" << std::endl;
251 std::clog <<
"Dimension : " << dimension << std::endl;
252 std::clog <<
"Euler characteristic : " << euler << std::endl;
253 std::clog <<
"Num simplices : " << num_simplices << std::endl;
254 std::clog <<
"Total time: " << clock << std::endl;
255 std::clog <<
"Time per simplex: " << clock.num_seconds() / num_simplices <<
" s" << std::endl;
256 std::clog <<
"+++++++++++++++++++++++++" << std::endl;
259 void show_homology_group() {
261 std::clog <<
"Works only on linux x64 for the moment\n";
269 void show_euler_characteristic() {
270 unsigned num_simplices = 0;
275 dimension = (std::max)(s.dimension(), dimension);
276 if (s.dimension() % 2 == 0)
281 std::clog <<
"Saw " << num_simplices <<
" simplices with maximum dimension " << dimension << std::endl;
282 std::clog <<
"The euler characteristic is : " << euler << std::endl;
285 void show_persistence(
int p,
double threshold,
int max_dim,
double min_pers) {
289 void show_critical_points(
double max_distance) {
293 void show_is_manifold() {
299 std::clog <<
"The complex is a " << dim <<
"-manifold\n";
302 std::clog <<
"The complex has dimension greater than " << dim <<
" and is not a manifold\n";
304 std::clog <<
"The complex has dimension>=4 and may or may not be a manifold\n";
311 save_complex_in_file_for_chomp();
312 std::clog <<
"Call CHOMP library\n";
313 int returnValue = system(
"homsimpl chomp.sim");
314 if (returnValue != 0) {
315 std::cerr <<
"homsimpl (from CHOMP) failed. Please check it is installed or available in the PATH."
320 void save_complex_in_file_for_chomp() {
322 file.open(
"chomp.sim");
338 unsigned num_vertices()
const {
339 return complex_.num_vertices();
342 unsigned num_edges()
const {
343 return complex_.num_edges();
Definition: Critical_points.h:26
Definition: Edge_collapsor.h:21
Definition: Edge_contractor.h:24
Edge_handle add_edge_without_blockers(Vertex_handle a, Vertex_handle b)
Adds an edge between vertices a and b without blockers.
Definition: Skeleton_blocker_complex.h:566
void keep_only_vertices()
The complex is reduced to its set of vertices. All the edges and blockers are removed.
Definition: Skeleton_blocker_complex.h:623
Complex_edge_range edge_range() const
Returns a Complex_edge_range over all edges of the simplicial complex.
Definition: Skeleton_blocker_complex.h:1318
Vertex_handle first_vertex(Edge_handle edge_handle) const
returns the first vertex of an edge
Definition: Skeleton_blocker_complex.h:512
int degree(Vertex_handle local) const
return the graph degree of a vertex.
Definition: Skeleton_blocker_complex.h:468
Vertex_handle second_vertex(Edge_handle edge_handle) const
returns the first vertex of an edge
Definition: Skeleton_blocker_complex.h:520
int num_connected_components() const
returns the number of connected components in the graph of the 1-skeleton.
Definition: Skeleton_blocker_complex.h:1029
virtual void clear()
Definition: Skeleton_blocker_complex.h:311
Complex_vertex_range vertex_range() const
Returns a Complex_vertex_range over all vertices of the complex.
Definition: Skeleton_blocker_complex.h:1284
Complex_simplex_range complex_simplex_range() const
Returns a Complex_simplex_range over all the simplices of the complex.
Definition: Skeleton_blocker_complex.h:1431
Vertex_handle add_vertex()
Add a vertex to the complex with a default constructed associated point.
Definition: Skeleton_blocker_geometric_complex.h:85
const Point & point(Vertex_handle v) const
Returns the Point associated to the vertex v.
Definition: Skeleton_blocker_geometric_complex.h:101
Definition: Is_manifold.h:24
Definition: Lloyd_builder.h:20
Definition: Persistence_compute.h:36
Definition: Vertex_collapsor.h:23
Definition: SkeletonBlockerDS.h:56
Handle type for the vertices of a cell complex.
Definition: VertexHandle.h:15