random_point_generators.h
1 /* This file is part of the Gudhi Library. The Gudhi library
2  * (Geometric Understanding in Higher Dimensions) is a generic C++
3  * library for computational topology.
4  *
5  * Author(s): Clement Jamin
6  *
7  * Copyright (C) 2016 Inria
8  *
9  * This program is free software: you can redistribute it and/or modify
10  * it under the terms of the GNU General Public License as published by
11  * the Free Software Foundation, either version 3 of the License, or
12  * (at your option) any later version.
13  *
14  * This program is distributed in the hope that it will be useful,
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17  * GNU General Public License for more details.
18  *
19  * You should have received a copy of the GNU General Public License
20  * along with this program. If not, see <http://www.gnu.org/licenses/>.
21  */
22 
23 #ifndef RANDOM_POINT_GENERATORS_H_
24 #define RANDOM_POINT_GENERATORS_H_
25 
26 #include <CGAL/number_utils.h>
27 #include <CGAL/Random.h>
28 #include <CGAL/point_generators_d.h>
29 
30 #include <vector> // for vector<>
31 
32 namespace Gudhi {
33 
35 // Note: All these functions have been tested with the CGAL::Epick_d kernel
37 
38 // construct_point: dim 2
39 
40 template <typename Kernel>
41 typename Kernel::Point_d construct_point(const Kernel &k,
42  typename Kernel::FT x1, typename Kernel::FT x2) {
43  typename Kernel::FT tab[2];
44  tab[0] = x1;
45  tab[1] = x2;
46  return k.construct_point_d_object()(2, &tab[0], &tab[2]);
47 }
48 
49 // construct_point: dim 3
50 
51 template <typename Kernel>
52 typename Kernel::Point_d construct_point(const Kernel &k,
53  typename Kernel::FT x1, typename Kernel::FT x2, typename Kernel::FT x3) {
54  typename Kernel::FT tab[3];
55  tab[0] = x1;
56  tab[1] = x2;
57  tab[2] = x3;
58  return k.construct_point_d_object()(3, &tab[0], &tab[3]);
59 }
60 
61 // construct_point: dim 4
62 
63 template <typename Kernel>
64 typename Kernel::Point_d construct_point(const Kernel &k,
65  typename Kernel::FT x1, typename Kernel::FT x2, typename Kernel::FT x3,
66  typename Kernel::FT x4) {
67  typename Kernel::FT tab[4];
68  tab[0] = x1;
69  tab[1] = x2;
70  tab[2] = x3;
71  tab[3] = x4;
72  return k.construct_point_d_object()(4, &tab[0], &tab[4]);
73 }
74 
75 // construct_point: dim 5
76 
77 template <typename Kernel>
78 typename Kernel::Point_d construct_point(const Kernel &k,
79  typename Kernel::FT x1, typename Kernel::FT x2, typename Kernel::FT x3,
80  typename Kernel::FT x4, typename Kernel::FT x5) {
81  typename Kernel::FT tab[5];
82  tab[0] = x1;
83  tab[1] = x2;
84  tab[2] = x3;
85  tab[3] = x4;
86  tab[4] = x5;
87  return k.construct_point_d_object()(5, &tab[0], &tab[5]);
88 }
89 
90 // construct_point: dim 6
91 
92 template <typename Kernel>
93 typename Kernel::Point_d construct_point(const Kernel &k,
94  typename Kernel::FT x1, typename Kernel::FT x2, typename Kernel::FT x3,
95  typename Kernel::FT x4, typename Kernel::FT x5, typename Kernel::FT x6) {
96  typename Kernel::FT tab[6];
97  tab[0] = x1;
98  tab[1] = x2;
99  tab[2] = x3;
100  tab[3] = x4;
101  tab[4] = x5;
102  tab[5] = x6;
103  return k.construct_point_d_object()(6, &tab[0], &tab[6]);
104 }
105 
106 template <typename Kernel>
107 std::vector<typename Kernel::Point_d> generate_points_on_plane(std::size_t num_points, int intrinsic_dim,
108  int ambient_dim,
109  double coord_min = -5., double coord_max = 5.) {
110  typedef typename Kernel::Point_d Point;
111  typedef typename Kernel::FT FT;
112  Kernel k;
113  CGAL::Random rng;
114  std::vector<Point> points;
115  points.reserve(num_points);
116  for (std::size_t i = 0; i < num_points;) {
117  std::vector<FT> pt(ambient_dim, FT(0));
118  for (int j = 0; j < intrinsic_dim; ++j)
119  pt[j] = rng.get_double(coord_min, coord_max);
120 
121  Point p = k.construct_point_d_object()(ambient_dim, pt.begin(), pt.end());
122  points.push_back(p);
123  ++i;
124  }
125  return points;
126 }
127 
128 template <typename Kernel>
129 std::vector<typename Kernel::Point_d> generate_points_on_moment_curve(std::size_t num_points, int dim,
130  typename Kernel::FT min_x,
131  typename Kernel::FT max_x) {
132  typedef typename Kernel::Point_d Point;
133  typedef typename Kernel::FT FT;
134  Kernel k;
135  CGAL::Random rng;
136  std::vector<Point> points;
137  points.reserve(num_points);
138  for (std::size_t i = 0; i < num_points;) {
139  FT x = rng.get_double(min_x, max_x);
140  std::vector<FT> coords;
141  coords.reserve(dim);
142  for (int p = 1; p <= dim; ++p)
143  coords.push_back(std::pow(CGAL::to_double(x), p));
144  Point p = k.construct_point_d_object()(
145  dim, coords.begin(), coords.end());
146  points.push_back(p);
147  ++i;
148  }
149  return points;
150 }
151 
152 
153 // R = big radius, r = small radius
154 template <typename Kernel/*, typename TC_basis*/>
155 std::vector<typename Kernel::Point_d> generate_points_on_torus_3D(std::size_t num_points, double R, double r,
156  bool uniform = false) {
157  typedef typename Kernel::Point_d Point;
158  typedef typename Kernel::FT FT;
159  Kernel k;
160  CGAL::Random rng;
161 
162  // if uniform
163  std::size_t num_lines = (std::size_t)sqrt(num_points);
164 
165  std::vector<Point> points;
166  points.reserve(num_points);
167  for (std::size_t i = 0; i < num_points;) {
168  FT u, v;
169  if (uniform) {
170  std::size_t k1 = i / num_lines;
171  std::size_t k2 = i % num_lines;
172  u = 6.2832 * k1 / num_lines;
173  v = 6.2832 * k2 / num_lines;
174  } else {
175  u = rng.get_double(0, 6.2832);
176  v = rng.get_double(0, 6.2832);
177  }
178  Point p = construct_point(k,
179  (R + r * std::cos(u)) * std::cos(v),
180  (R + r * std::cos(u)) * std::sin(v),
181  r * std::sin(u));
182  points.push_back(p);
183  ++i;
184  }
185  return points;
186 }
187 
188 // "Private" function used by generate_points_on_torus_d
189 template <typename Kernel, typename OutputIterator>
190 static void generate_uniform_points_on_torus_d(const Kernel &k, int dim, std::size_t num_slices,
191  OutputIterator out,
192  double radius_noise_percentage = 0.,
193  std::vector<typename Kernel::FT> current_point =
194  std::vector<typename Kernel::FT>()) {
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 = 6.2832 * 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_uniform_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, bool uniform = false,
219  double radius_noise_percentage = 0.) {
220  typedef typename Kernel::Point_d Point;
221  typedef typename Kernel::FT FT;
222  Kernel k;
223  CGAL::Random rng;
224 
225  std::vector<Point> points;
226  points.reserve(num_points);
227  if (uniform) {
228  std::size_t num_slices = (std::size_t)std::pow(num_points, 1. / dim);
229  generate_uniform_points_on_torus_d(
230  k, dim, num_slices, std::back_inserter(points), radius_noise_percentage);
231  } else {
232  for (std::size_t i = 0; i < num_points;) {
233  double radius_noise_ratio = 1.;
234  if (radius_noise_percentage > 0.) {
235  radius_noise_ratio = rng.get_double(
236  (100. - radius_noise_percentage) / 100.,
237  (100. + radius_noise_percentage) / 100.);
238  }
239  std::vector<typename Kernel::FT> pt;
240  pt.reserve(dim * 2);
241  for (int curdim = 0; curdim < dim; ++curdim) {
242  FT alpha = rng.get_double(0, 6.2832);
243  pt.push_back(radius_noise_ratio * std::cos(alpha));
244  pt.push_back(radius_noise_ratio * std::sin(alpha));
245  }
246 
247  Point p = k.construct_point_d_object()(pt.begin(), pt.end());
248  points.push_back(p);
249  ++i;
250  }
251  }
252  return points;
253 }
254 
255 template <typename Kernel>
256 std::vector<typename Kernel::Point_d> generate_points_on_sphere_d(std::size_t num_points, int dim, double radius,
257  double radius_noise_percentage = 0.) {
258  typedef typename Kernel::Point_d Point;
259  Kernel k;
260  CGAL::Random rng;
261  CGAL::Random_points_on_sphere_d<Point> generator(dim, radius);
262  std::vector<Point> points;
263  points.reserve(num_points);
264  for (std::size_t i = 0; i < num_points;) {
265  Point p = *generator++;
266  if (radius_noise_percentage > 0.) {
267  double radius_noise_ratio = rng.get_double(
268  (100. - radius_noise_percentage) / 100.,
269  (100. + radius_noise_percentage) / 100.);
270 
271  typename Kernel::Point_to_vector_d k_pt_to_vec =
272  k.point_to_vector_d_object();
273  typename Kernel::Vector_to_point_d k_vec_to_pt =
274  k.vector_to_point_d_object();
275  typename Kernel::Scaled_vector_d k_scaled_vec =
276  k.scaled_vector_d_object();
277  p = k_vec_to_pt(k_scaled_vec(k_pt_to_vec(p), radius_noise_ratio));
278  }
279  points.push_back(p);
280  ++i;
281  }
282  return points;
283 }
284 
285 template <typename Kernel>
286 std::vector<typename Kernel::Point_d> generate_points_in_ball_d(std::size_t num_points, int dim, double radius) {
287  typedef typename Kernel::Point_d Point;
288  Kernel k;
289  CGAL::Random rng;
290  CGAL::Random_points_in_ball_d<Point> generator(dim, radius);
291  std::vector<Point> points;
292  points.reserve(num_points);
293  for (std::size_t i = 0; i < num_points;) {
294  Point p = *generator++;
295  points.push_back(p);
296  ++i;
297  }
298  return points;
299 }
300 
301 template <typename Kernel>
302 std::vector<typename Kernel::Point_d> generate_points_in_cube_d(std::size_t num_points, int dim, double radius) {
303  typedef typename Kernel::Point_d Point;
304  Kernel k;
305  CGAL::Random rng;
306  CGAL::Random_points_in_cube_d<Point> generator(dim, radius);
307  std::vector<Point> points;
308  points.reserve(num_points);
309  for (std::size_t i = 0; i < num_points;) {
310  Point p = *generator++;
311  points.push_back(p);
312  ++i;
313  }
314  return points;
315 }
316 
317 template <typename Kernel>
318 std::vector<typename Kernel::Point_d> generate_points_on_two_spheres_d(std::size_t num_points, int dim, double radius,
319  double distance_between_centers,
320  double radius_noise_percentage = 0.) {
321  typedef typename Kernel::FT FT;
322  typedef typename Kernel::Point_d Point;
323  typedef typename Kernel::Vector_d Vector;
324  Kernel k;
325  CGAL::Random rng;
326  CGAL::Random_points_on_sphere_d<Point> generator(dim, radius);
327  std::vector<Point> points;
328  points.reserve(num_points);
329 
330  std::vector<FT> t(dim, FT(0));
331  t[0] = distance_between_centers;
332  Vector c1_to_c2(t.begin(), t.end());
333 
334  for (std::size_t i = 0; i < num_points;) {
335  Point p = *generator++;
336  if (radius_noise_percentage > 0.) {
337  double radius_noise_ratio = rng.get_double(
338  (100. - radius_noise_percentage) / 100.,
339  (100. + radius_noise_percentage) / 100.);
340 
341  typename Kernel::Point_to_vector_d k_pt_to_vec =
342  k.point_to_vector_d_object();
343  typename Kernel::Vector_to_point_d k_vec_to_pt =
344  k.vector_to_point_d_object();
345  typename Kernel::Scaled_vector_d k_scaled_vec =
346  k.scaled_vector_d_object();
347  p = k_vec_to_pt(k_scaled_vec(k_pt_to_vec(p), radius_noise_ratio));
348  }
349 
350  typename Kernel::Translated_point_d k_transl =
351  k.translated_point_d_object();
352  Point p2 = k_transl(p, c1_to_c2);
353  points.push_back(p);
354  points.push_back(p2);
355  i += 2;
356  }
357  return points;
358 }
359 
360 // Product of a 3-sphere and a circle => d = 3 / D = 5
361 
362 template <typename Kernel>
363 std::vector<typename Kernel::Point_d> generate_points_on_3sphere_and_circle(std::size_t num_points,
364  double sphere_radius) {
365  typedef typename Kernel::FT FT;
366  typedef typename Kernel::Point_d Point;
367  Kernel k;
368  CGAL::Random rng;
369  CGAL::Random_points_on_sphere_d<Point> generator(3, sphere_radius);
370  std::vector<Point> points;
371  points.reserve(num_points);
372 
373  typename Kernel::Compute_coordinate_d k_coord =
374  k.compute_coordinate_d_object();
375  for (std::size_t i = 0; i < num_points;) {
376  Point p_sphere = *generator++; // First 3 coords
377 
378  FT alpha = rng.get_double(0, 6.2832);
379  std::vector<FT> pt(5);
380  pt[0] = k_coord(p_sphere, 0);
381  pt[1] = k_coord(p_sphere, 1);
382  pt[2] = k_coord(p_sphere, 2);
383  pt[3] = std::cos(alpha);
384  pt[4] = std::sin(alpha);
385  Point p(pt.begin(), pt.end());
386  points.push_back(p);
387  ++i;
388  }
389  return points;
390 }
391 
392 // a = big radius, b = small radius
393 template <typename Kernel>
394 std::vector<typename Kernel::Point_d> generate_points_on_klein_bottle_3D(std::size_t num_points, double a, double b,
395  bool uniform = false) {
396  typedef typename Kernel::Point_d Point;
397  typedef typename Kernel::FT FT;
398  Kernel k;
399  CGAL::Random rng;
400 
401  // if uniform
402  std::size_t num_lines = (std::size_t)sqrt(num_points);
403 
404  std::vector<Point> points;
405  points.reserve(num_points);
406  for (std::size_t i = 0; i < num_points;) {
407  FT u, v;
408  if (uniform) {
409  std::size_t k1 = i / num_lines;
410  std::size_t k2 = i % num_lines;
411  u = 6.2832 * k1 / num_lines;
412  v = 6.2832 * k2 / num_lines;
413  } else {
414  u = rng.get_double(0, 6.2832);
415  v = rng.get_double(0, 6.2832);
416  }
417  double tmp = cos(u / 2) * sin(v) - sin(u / 2) * sin(2. * v);
418  Point p = construct_point(k,
419  (a + b * tmp) * cos(u),
420  (a + b * tmp) * sin(u),
421  b * (sin(u / 2) * sin(v) + cos(u / 2) * sin(2. * v)));
422  points.push_back(p);
423  ++i;
424  }
425  return points;
426 }
427 
428 // a = big radius, b = small radius
429 template <typename Kernel>
430 std::vector<typename Kernel::Point_d> generate_points_on_klein_bottle_4D(std::size_t num_points, double a, double b,
431  double noise = 0., bool uniform = false) {
432  typedef typename Kernel::Point_d Point;
433  typedef typename Kernel::FT FT;
434  Kernel k;
435  CGAL::Random rng;
436 
437  // if uniform
438  std::size_t num_lines = (std::size_t)sqrt(num_points);
439 
440  std::vector<Point> points;
441  points.reserve(num_points);
442  for (std::size_t i = 0; i < num_points;) {
443  FT u, v;
444  if (uniform) {
445  std::size_t k1 = i / num_lines;
446  std::size_t k2 = i % num_lines;
447  u = 6.2832 * k1 / num_lines;
448  v = 6.2832 * k2 / num_lines;
449  } else {
450  u = rng.get_double(0, 6.2832);
451  v = rng.get_double(0, 6.2832);
452  }
453  Point p = construct_point(k,
454  (a + b * cos(v)) * cos(u) + (noise == 0. ? 0. : rng.get_double(0, noise)),
455  (a + b * cos(v)) * sin(u) + (noise == 0. ? 0. : rng.get_double(0, noise)),
456  b * sin(v) * cos(u / 2) + (noise == 0. ? 0. : rng.get_double(0, noise)),
457  b * sin(v) * sin(u / 2) + (noise == 0. ? 0. : rng.get_double(0, noise)));
458  points.push_back(p);
459  ++i;
460  }
461  return points;
462 }
463 
464 
465 // a = big radius, b = small radius
466 
467 template <typename Kernel>
468 std::vector<typename Kernel::Point_d>
469 generate_points_on_klein_bottle_variant_5D(
470  std::size_t num_points, double a, double b, bool uniform = false) {
471  typedef typename Kernel::Point_d Point;
472  typedef typename Kernel::FT FT;
473  Kernel k;
474  CGAL::Random rng;
475 
476  // if uniform
477  std::size_t num_lines = (std::size_t)sqrt(num_points);
478 
479  std::vector<Point> points;
480  points.reserve(num_points);
481  for (std::size_t i = 0; i < num_points;) {
482  FT u, v;
483  if (uniform) {
484  std::size_t k1 = i / num_lines;
485  std::size_t k2 = i % num_lines;
486  u = 6.2832 * k1 / num_lines;
487  v = 6.2832 * k2 / num_lines;
488  } else {
489  u = rng.get_double(0, 6.2832);
490  v = rng.get_double(0, 6.2832);
491  }
492  FT x1 = (a + b * cos(v)) * cos(u);
493  FT x2 = (a + b * cos(v)) * sin(u);
494  FT x3 = b * sin(v) * cos(u / 2);
495  FT x4 = b * sin(v) * sin(u / 2);
496  FT x5 = x1 + x2 + x3 + x4;
497 
498  Point p = construct_point(k, x1, x2, x3, x4, x5);
499  points.push_back(p);
500  ++i;
501  }
502  return points;
503 }
504 
505 } // namespace Gudhi
506 
507 #endif // RANDOM_POINT_GENERATORS_H_
Definition: SimplicialComplexForAlpha.h:26
GUDHI  Version 2.2.0  - C++ library for Topological Data Analysis (TDA) and Higher Dimensional Geometry Understanding.  - Copyright : GPL v3 Generated on Tue Aug 14 2018 11:45:43 for GUDHI by Doxygen 1.8.13