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