random_point_generators.h
1 /* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
2  * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
3  * Author(s): Clement Jamin
4  *
5  * Copyright (C) 2016 Inria
6  *
7  * Modification(s):
8  * - 2019/08 Vincent Rouvreau: Fix issue #10 for CGAL
9  * - YYYY/MM Author: Description of the modification
10  */
11 
12 #ifndef RANDOM_POINT_GENERATORS_H_
13 #define RANDOM_POINT_GENERATORS_H_
14 
15 #include <CGAL/number_utils.h>
16 #include <CGAL/Random.h>
17 #include <CGAL/point_generators_d.h>
18 #include <CGAL/version.h> // for CGAL_VERSION_NR
19 
20 #include <vector> // for vector<>
21 #include <boost/math/constants/constants.hpp> // for pi constant
22 
23 // Make compilation fail - required for external projects - https://github.com/GUDHI/gudhi-devel/issues/10
24 #if CGAL_VERSION_NR < 1041101000
25 # error random_point_generators is only available for CGAL >= 4.11
26 #endif
27 
28 namespace Gudhi {
29 
31 // Note: All these functions have been tested with the CGAL::Epick_d kernel
33 
34 // construct_point: dim 2
35 
36 template <typename Kernel>
37 typename Kernel::Point_d construct_point(const Kernel &k,
38  typename Kernel::FT x1, typename Kernel::FT x2) {
39  typename Kernel::FT tab[2];
40  tab[0] = x1;
41  tab[1] = x2;
42  return k.construct_point_d_object()(2, &tab[0], &tab[2]);
43 }
44 
45 // construct_point: dim 3
46 
47 template <typename Kernel>
48 typename Kernel::Point_d construct_point(const Kernel &k,
49  typename Kernel::FT x1, typename Kernel::FT x2, typename Kernel::FT x3) {
50  typename Kernel::FT tab[3];
51  tab[0] = x1;
52  tab[1] = x2;
53  tab[2] = x3;
54  return k.construct_point_d_object()(3, &tab[0], &tab[3]);
55 }
56 
57 // construct_point: dim 4
58 
59 template <typename Kernel>
60 typename Kernel::Point_d construct_point(const Kernel &k,
61  typename Kernel::FT x1, typename Kernel::FT x2, typename Kernel::FT x3,
62  typename Kernel::FT x4) {
63  typename Kernel::FT tab[4];
64  tab[0] = x1;
65  tab[1] = x2;
66  tab[2] = x3;
67  tab[3] = x4;
68  return k.construct_point_d_object()(4, &tab[0], &tab[4]);
69 }
70 
71 // construct_point: dim 5
72 
73 template <typename Kernel>
74 typename Kernel::Point_d construct_point(const Kernel &k,
75  typename Kernel::FT x1, typename Kernel::FT x2, typename Kernel::FT x3,
76  typename Kernel::FT x4, typename Kernel::FT x5) {
77  typename Kernel::FT tab[5];
78  tab[0] = x1;
79  tab[1] = x2;
80  tab[2] = x3;
81  tab[3] = x4;
82  tab[4] = x5;
83  return k.construct_point_d_object()(5, &tab[0], &tab[5]);
84 }
85 
86 // construct_point: dim 6
87 
88 template <typename Kernel>
89 typename Kernel::Point_d construct_point(const Kernel &k,
90  typename Kernel::FT x1, typename Kernel::FT x2, typename Kernel::FT x3,
91  typename Kernel::FT x4, typename Kernel::FT x5, typename Kernel::FT x6) {
92  typename Kernel::FT tab[6];
93  tab[0] = x1;
94  tab[1] = x2;
95  tab[2] = x3;
96  tab[3] = x4;
97  tab[4] = x5;
98  tab[5] = x6;
99  return k.construct_point_d_object()(6, &tab[0], &tab[6]);
100 }
101 
102 template <typename Kernel>
103 std::vector<typename Kernel::Point_d> generate_points_on_plane(std::size_t num_points, int intrinsic_dim,
104  int ambient_dim,
105  double coord_min = -5., double coord_max = 5.) {
106  typedef typename Kernel::Point_d Point;
107  typedef typename Kernel::FT FT;
108  Kernel k;
109  CGAL::Random rng;
110  std::vector<Point> points;
111  points.reserve(num_points);
112  for (std::size_t i = 0; i < num_points;) {
113  std::vector<FT> pt(ambient_dim, FT(0));
114  for (int j = 0; j < intrinsic_dim; ++j)
115  pt[j] = rng.get_double(coord_min, coord_max);
116 
117  Point p = k.construct_point_d_object()(ambient_dim, pt.begin(), pt.end());
118  points.push_back(p);
119  ++i;
120  }
121  return points;
122 }
123 
124 template <typename Kernel>
125 std::vector<typename Kernel::Point_d> generate_points_on_moment_curve(std::size_t num_points, int dim,
126  typename Kernel::FT min_x,
127  typename Kernel::FT max_x) {
128  typedef typename Kernel::Point_d Point;
129  typedef typename Kernel::FT FT;
130  Kernel k;
131  CGAL::Random rng;
132  std::vector<Point> points;
133  points.reserve(num_points);
134  for (std::size_t i = 0; i < num_points;) {
135  FT x = rng.get_double(min_x, max_x);
136  std::vector<FT> coords;
137  coords.reserve(dim);
138  for (int p = 1; p <= dim; ++p)
139  coords.push_back(std::pow(CGAL::to_double(x), p));
140  Point p = k.construct_point_d_object()(
141  dim, coords.begin(), coords.end());
142  points.push_back(p);
143  ++i;
144  }
145  return points;
146 }
147 
148 
149 // R = big radius, r = small radius
150 template <typename Kernel/*, typename TC_basis*/>
151 std::vector<typename Kernel::Point_d> generate_points_on_torus_3D(std::size_t num_points, double R, double r,
152  bool uniform = false) {
153  using namespace boost::math::double_constants;
154 
155  typedef typename Kernel::Point_d Point;
156  typedef typename Kernel::FT FT;
157  Kernel k;
158  CGAL::Random rng;
159 
160  // if uniform
161  std::size_t num_lines = (std::size_t)sqrt(num_points);
162 
163  std::vector<Point> points;
164  points.reserve(num_points);
165  for (std::size_t i = 0; i < num_points;) {
166  FT u, v;
167  if (uniform) {
168  std::size_t k1 = i / num_lines;
169  std::size_t k2 = i % num_lines;
170  u = two_pi * k1 / num_lines;
171  v = two_pi * k2 / num_lines;
172  } else {
173  u = rng.get_double(0, two_pi);
174  v = rng.get_double(0, two_pi);
175  }
176  Point p = construct_point(k,
177  (R + r * std::cos(u)) * std::cos(v),
178  (R + r * std::cos(u)) * std::sin(v),
179  r * std::sin(u));
180  points.push_back(p);
181  ++i;
182  }
183  return points;
184 }
185 
186 // "Private" function used by generate_points_on_torus_d
187 template <typename Kernel, typename OutputIterator>
188 static void generate_grid_points_on_torus_d(const Kernel &k, int dim, std::size_t num_slices,
189  OutputIterator out,
190  double radius_noise_percentage = 0.,
191  std::vector<typename Kernel::FT> current_point =
192  std::vector<typename Kernel::FT>()) {
193  using namespace boost::math::double_constants;
194 
195  CGAL::Random rng;
196  int point_size = static_cast<int>(current_point.size());
197  if (point_size == 2 * dim) {
198  *out++ = k.construct_point_d_object()(point_size, current_point.begin(), current_point.end());
199  } else {
200  for (std::size_t slice_idx = 0; slice_idx < num_slices; ++slice_idx) {
201  double radius_noise_ratio = 1.;
202  if (radius_noise_percentage > 0.) {
203  radius_noise_ratio = rng.get_double(
204  (100. - radius_noise_percentage) / 100.,
205  (100. + radius_noise_percentage) / 100.);
206  }
207  std::vector<typename Kernel::FT> cp2 = current_point;
208  double alpha = two_pi * slice_idx / num_slices;
209  cp2.push_back(radius_noise_ratio * std::cos(alpha));
210  cp2.push_back(radius_noise_ratio * std::sin(alpha));
211  generate_grid_points_on_torus_d(
212  k, dim, num_slices, out, radius_noise_percentage, cp2);
213  }
214  }
215 }
216 
217 template <typename Kernel>
218 std::vector<typename Kernel::Point_d> generate_points_on_torus_d(std::size_t num_points, int dim, std::string sample = "random",
219  double radius_noise_percentage = 0.) {
220  using namespace boost::math::double_constants;
221 
222  typedef typename Kernel::Point_d Point;
223  typedef typename Kernel::FT FT;
224  Kernel k;
225  CGAL::Random rng;
226 
227  std::vector<Point> points;
228  points.reserve(num_points);
229  if (sample == "grid") {
230  std::size_t num_slices = (std::size_t)std::pow(num_points + .5, 1. / dim); // add .5 to avoid rounding down with numerical approximations
231  generate_grid_points_on_torus_d(
232  k, dim, num_slices, std::back_inserter(points), radius_noise_percentage);
233  } else {
234  for (std::size_t i = 0; i < num_points;) {
235  double radius_noise_ratio = 1.;
236  if (radius_noise_percentage > 0.) {
237  radius_noise_ratio = rng.get_double(
238  (100. - radius_noise_percentage) / 100.,
239  (100. + radius_noise_percentage) / 100.);
240  }
241  std::vector<typename Kernel::FT> pt;
242  pt.reserve(dim * 2);
243  for (int curdim = 0; curdim < dim; ++curdim) {
244  FT alpha = rng.get_double(0, two_pi);
245  pt.push_back(radius_noise_ratio * std::cos(alpha));
246  pt.push_back(radius_noise_ratio * std::sin(alpha));
247  }
248 
249  Point p = k.construct_point_d_object()(pt.begin(), pt.end());
250  points.push_back(p);
251  ++i;
252  }
253  }
254  return points;
255 }
256 
257 template <typename Kernel>
258 std::vector<typename Kernel::Point_d> generate_points_on_sphere_d(std::size_t num_points, int dim, double radius,
259  double radius_noise_percentage = 0.) {
260  typedef typename Kernel::Point_d Point;
261  Kernel k;
262  CGAL::Random rng;
263  CGAL::Random_points_on_sphere_d<Point> generator(dim, radius);
264  std::vector<Point> points;
265  points.reserve(num_points);
266  for (std::size_t i = 0; i < num_points;) {
267  Point p = *generator++;
268  if (radius_noise_percentage > 0.) {
269  double radius_noise_ratio = rng.get_double(
270  (100. - radius_noise_percentage) / 100.,
271  (100. + radius_noise_percentage) / 100.);
272 
273  typename Kernel::Point_to_vector_d k_pt_to_vec =
274  k.point_to_vector_d_object();
275  typename Kernel::Vector_to_point_d k_vec_to_pt =
276  k.vector_to_point_d_object();
277  typename Kernel::Scaled_vector_d k_scaled_vec =
278  k.scaled_vector_d_object();
279  p = k_vec_to_pt(k_scaled_vec(k_pt_to_vec(p), radius_noise_ratio));
280  }
281  points.push_back(p);
282  ++i;
283  }
284  return points;
285 }
286 
287 template <typename Kernel>
288 std::vector<typename Kernel::Point_d> generate_points_in_ball_d(std::size_t num_points, int dim, double radius) {
289  typedef typename Kernel::Point_d Point;
290  Kernel k;
291  CGAL::Random rng;
292  CGAL::Random_points_in_ball_d<Point> generator(dim, radius);
293  std::vector<Point> points;
294  points.reserve(num_points);
295  for (std::size_t i = 0; i < num_points;) {
296  Point p = *generator++;
297  points.push_back(p);
298  ++i;
299  }
300  return points;
301 }
302 
303 template <typename Kernel>
304 std::vector<typename Kernel::Point_d> generate_points_in_cube_d(std::size_t num_points, int dim, double radius) {
305  typedef typename Kernel::Point_d Point;
306  Kernel k;
307  CGAL::Random rng;
308  CGAL::Random_points_in_cube_d<Point> generator(dim, radius);
309  std::vector<Point> points;
310  points.reserve(num_points);
311  for (std::size_t i = 0; i < num_points;) {
312  Point p = *generator++;
313  points.push_back(p);
314  ++i;
315  }
316  return points;
317 }
318 
319 template <typename Kernel>
320 std::vector<typename Kernel::Point_d> generate_points_on_two_spheres_d(std::size_t num_points, int dim, double radius,
321  double distance_between_centers,
322  double radius_noise_percentage = 0.) {
323  typedef typename Kernel::FT FT;
324  typedef typename Kernel::Point_d Point;
325  typedef typename Kernel::Vector_d Vector;
326  Kernel k;
327  CGAL::Random rng;
328  CGAL::Random_points_on_sphere_d<Point> generator(dim, radius);
329  std::vector<Point> points;
330  points.reserve(num_points);
331 
332  std::vector<FT> t(dim, FT(0));
333  t[0] = distance_between_centers;
334  Vector c1_to_c2(t.begin(), t.end());
335 
336  for (std::size_t i = 0; i < num_points;) {
337  Point p = *generator++;
338  if (radius_noise_percentage > 0.) {
339  double radius_noise_ratio = rng.get_double(
340  (100. - radius_noise_percentage) / 100.,
341  (100. + radius_noise_percentage) / 100.);
342 
343  typename Kernel::Point_to_vector_d k_pt_to_vec =
344  k.point_to_vector_d_object();
345  typename Kernel::Vector_to_point_d k_vec_to_pt =
346  k.vector_to_point_d_object();
347  typename Kernel::Scaled_vector_d k_scaled_vec =
348  k.scaled_vector_d_object();
349  p = k_vec_to_pt(k_scaled_vec(k_pt_to_vec(p), radius_noise_ratio));
350  }
351 
352  typename Kernel::Translated_point_d k_transl =
353  k.translated_point_d_object();
354  Point p2 = k_transl(p, c1_to_c2);
355  points.push_back(p);
356  points.push_back(p2);
357  i += 2;
358  }
359  return points;
360 }
361 
362 // Product of a 3-sphere and a circle => d = 3 / D = 5
363 
364 template <typename Kernel>
365 std::vector<typename Kernel::Point_d> generate_points_on_3sphere_and_circle(std::size_t num_points,
366  double sphere_radius) {
367  using namespace boost::math::double_constants;
368 
369  typedef typename Kernel::FT FT;
370  typedef typename Kernel::Point_d Point;
371  Kernel k;
372  CGAL::Random rng;
373  CGAL::Random_points_on_sphere_d<Point> generator(3, sphere_radius);
374  std::vector<Point> points;
375  points.reserve(num_points);
376 
377  typename Kernel::Compute_coordinate_d k_coord =
378  k.compute_coordinate_d_object();
379  for (std::size_t i = 0; i < num_points;) {
380  Point p_sphere = *generator++; // First 3 coords
381 
382  FT alpha = rng.get_double(0, two_pi);
383  std::vector<FT> pt(5);
384  pt[0] = k_coord(p_sphere, 0);
385  pt[1] = k_coord(p_sphere, 1);
386  pt[2] = k_coord(p_sphere, 2);
387  pt[3] = std::cos(alpha);
388  pt[4] = std::sin(alpha);
389  Point p(pt.begin(), pt.end());
390  points.push_back(p);
391  ++i;
392  }
393  return points;
394 }
395 
396 // a = big radius, b = small radius
397 template <typename Kernel>
398 std::vector<typename Kernel::Point_d> generate_points_on_klein_bottle_3D(std::size_t num_points, double a, double b,
399  bool uniform = false) {
400  using namespace boost::math::double_constants;
401 
402  typedef typename Kernel::Point_d Point;
403  typedef typename Kernel::FT FT;
404  Kernel k;
405  CGAL::Random rng;
406 
407  // if uniform
408  std::size_t num_lines = (std::size_t)sqrt(num_points);
409 
410  std::vector<Point> points;
411  points.reserve(num_points);
412  for (std::size_t i = 0; i < num_points;) {
413  FT u, v;
414  if (uniform) {
415  std::size_t k1 = i / num_lines;
416  std::size_t k2 = i % num_lines;
417  u = two_pi * k1 / num_lines;
418  v = two_pi * k2 / num_lines;
419  } else {
420  u = rng.get_double(0, two_pi);
421  v = rng.get_double(0, two_pi);
422  }
423  double tmp = cos(u / 2) * sin(v) - sin(u / 2) * sin(2. * v);
424  Point p = construct_point(k,
425  (a + b * tmp) * cos(u),
426  (a + b * tmp) * sin(u),
427  b * (sin(u / 2) * sin(v) + cos(u / 2) * sin(2. * v)));
428  points.push_back(p);
429  ++i;
430  }
431  return points;
432 }
433 
434 // a = big radius, b = small radius
435 template <typename Kernel>
436 std::vector<typename Kernel::Point_d> generate_points_on_klein_bottle_4D(std::size_t num_points, double a, double b,
437  double noise = 0., bool uniform = false) {
438  using namespace boost::math::double_constants;
439 
440  typedef typename Kernel::Point_d Point;
441  typedef typename Kernel::FT FT;
442  Kernel k;
443  CGAL::Random rng;
444 
445  // if uniform
446  std::size_t num_lines = (std::size_t)sqrt(num_points);
447 
448  std::vector<Point> points;
449  points.reserve(num_points);
450  for (std::size_t i = 0; i < num_points;) {
451  FT u, v;
452  if (uniform) {
453  std::size_t k1 = i / num_lines;
454  std::size_t k2 = i % num_lines;
455  u = two_pi * k1 / num_lines;
456  v = two_pi * k2 / num_lines;
457  } else {
458  u = rng.get_double(0, two_pi);
459  v = rng.get_double(0, two_pi);
460  }
461  Point p = construct_point(k,
462  (a + b * cos(v)) * cos(u) + (noise == 0. ? 0. : rng.get_double(0, noise)),
463  (a + b * cos(v)) * sin(u) + (noise == 0. ? 0. : rng.get_double(0, noise)),
464  b * sin(v) * cos(u / 2) + (noise == 0. ? 0. : rng.get_double(0, noise)),
465  b * sin(v) * sin(u / 2) + (noise == 0. ? 0. : rng.get_double(0, noise)));
466  points.push_back(p);
467  ++i;
468  }
469  return points;
470 }
471 
472 
473 // a = big radius, b = small radius
474 
475 template <typename Kernel>
476 std::vector<typename Kernel::Point_d>
477 generate_points_on_klein_bottle_variant_5D(
478  std::size_t num_points, double a, double b, bool uniform = false) {
479  using namespace boost::math::double_constants;
480 
481  typedef typename Kernel::Point_d Point;
482  typedef typename Kernel::FT FT;
483  Kernel k;
484  CGAL::Random rng;
485 
486  // if uniform
487  std::size_t num_lines = (std::size_t)sqrt(num_points);
488 
489  std::vector<Point> points;
490  points.reserve(num_points);
491  for (std::size_t i = 0; i < num_points;) {
492  FT u, v;
493  if (uniform) {
494  std::size_t k1 = i / num_lines;
495  std::size_t k2 = i % num_lines;
496  u = two_pi * k1 / num_lines;
497  v = two_pi * k2 / num_lines;
498  } else {
499  u = rng.get_double(0, two_pi);
500  v = rng.get_double(0, two_pi);
501  }
502  FT x1 = (a + b * cos(v)) * cos(u);
503  FT x2 = (a + b * cos(v)) * sin(u);
504  FT x3 = b * sin(v) * cos(u / 2);
505  FT x4 = b * sin(v) * sin(u / 2);
506  FT x5 = x1 + x2 + x3 + x4;
507 
508  Point p = construct_point(k, x1, x2, x3, x4, x5);
509  points.push_back(p);
510  ++i;
511  }
512  return points;
513 }
514 
515 } // namespace Gudhi
516 
517 #endif // RANDOM_POINT_GENERATORS_H_
GUDHIdev  Version 3.5.0  - C++ library for Topological Data Analysis (TDA) and Higher Dimensional Geometry Understanding.  - Copyright : MIT Generated on Wed Apr 6 2022 19:26:28 for GUDHIdev by Doxygen 1.9.1